Numpyだけで気温と降水量の回帰分析【Python機械学習】
今回もまたNumpyだけを使った(sklearnなどを使わない)回帰分析。
実際の株価気象のデータを分析する。
実行環境
Androidスマホ
termux
Python3.8
JupyterNotebook
- Pythonライブラリ
- Numpy
- Pandas
- matplotlib
目次
目的
今回の目的も、Numpyのnp.polyfit()
とnp.poly1d()
をとりあえず使うとどんなことが出来るのかを実際のデータを使って体験します。
何の道具か、何の為に使うのか、何が戻り値なのか等、分かってないままやっています。
とりあえず使う。
題材として株価の日毎気象のデータを使ってグラフ作成をします。
- X軸 : 月別
- Y軸 :
株価気温、または降水量
そのデータを、Numpyのnp.polyfit()
とnp.poly1d()
にゴニョゴニョすると、気温や降水量の変動をなぞるような形を示す方程式が良い感じで導出できるようなので、簡単に試みます。
CSVデータの入手と注意
追記
株価データを使おうと途中まで作業していたんですが、「著作権法で云々」のためデータの転載ができないので別の部門に切り替える。使えない!
グラフも駄目なんだろうか?
【データはこちら】
日経平均プロフィルダウンロードセンター
ダウンロードセンター - 日経平均プロフィルここの「日経平均株価 日次データ csv」をダウンロードした。
次のファイルがダウンロードされる。
"nikkei_stock_average_daily_jp.csv"
追記
上記サイトのデータの転記等は著作権法により禁止されているのでブログで使っては駄目らしい。
株価データの転載に関するブログ記事。これを読むと殆どの株価データは使えない
なぁと。
無難に気象庁のお天気データでも使うとしますか。
国土交通省気象庁 過去の気象データダウンロードページ
気象庁|過去の気象データ・ダウンロード
特に著作権が云々と書いてないようなので大丈夫だろう。国の公表データだしね。
上のページで、順番に
- 「地点を選ぶ」
- 「項目を選ぶ」
- 「期間を選ぶ」
- 「CSVファイルをダウンロード」
で選択項目にチェック入れて「ダウンロード」すると、データが1つのCSVファイル「data.csv」の名前でダウンロードされます。
ブラウザのCookieをONにしないとできません。
日毎のデータを何十年分などデータ量が多くなると警告がでます。そういう時は小分けで。
長寿日本一、二を競い、健康寿命も長い長野県。諏訪大社の御加護でしょうか、食べ物や空気、水のおかげかな?気候がよろしい?
長脛彦などの「長」は蛇を表すそうでパイソンなだけに縁がありそうな、長野県諏訪の気象データを使います。
使うデータ
- ジャンル : 気象データ
- 対象地区 : 長野県諏訪
- 期間 : 1970年元旦から2020年元旦
- 種類 : 月別平均気温/月別降水量
Numpy.genfromtxt()でCSVファイルを読み込む(☠️)
参考
NumPyでCSVファイルを読み込み・書き込み(入力・出力) | note.nkmk.me
Numpyを使ったCSVファイルの読み込みにはいくつか方法があるようですが、今回はnp.genfromtxt()
を使ってみます。
特徴としては
- 複雑なCSVデータの読み込み
- 欠損値の処理
- 異なるデータ型の処理
ができる。
import numpy as np # 気象データのファイル data_file = "data.csv" # NumpyでCSVファイル読み込み data =np.genfromtxt( data_file, #names=True, skip_header=1, delimiter=',', encoding="shift-jis", dtype=None, #skip_footer=1 ) data
array([['', '諏訪', '諏訪', ..., '諏訪', '諏訪', '諏訪'],
['年月', '平均気温(℃)', '平均気温(℃)', ..., '降水量の合計(mm)', '降水量の合計(mm)',
'降水量の合計(mm)'],
['', '', '品質情報', ..., '現象なし情報', '品質情報', '均質番号'],
...,
['2019/11', '7.7', '8', ..., '0', '8', '1'],
['2019/12', '3.3', '8', ..., '0', '8', '1'],
['2020/1', '1.4', '4', ..., '0', '4', '1']], dtype='<U10')
CSVの読み込みでエラー出しまくる。パラメータをいじりまわす。エラーが出ない設定をつかむのが難しい。
上のリンク先「何度も言うがPandasでやれよ、簡単だぞ?!」
とりあえず読めたようなので何か確認する。
# データの型
data.dtype
dtype('<U10')
# 配列の形状
data.shape
(604, 8)
# 配列の要素数
data.size
4832
その後も色々とエラー出まくりング。諦めた。
Pandasでやろう。
NumpyでCSV読込は諦めてPandas.read_csv()で!
pandas.read_csv(CSVファイル, encoding="ナントカ")
でCSVデータがPandasのデータフレーム形式に変換される。
まずは無心で読み込む。
import pandas as pd # 気象庁よりDLした長野県諏訪のデータ data_file = "data.csv" # Pandasでcsvファイルの読込みD.F化 df_test = pd.read_csv( data_file, encoding = "shiftjis", ) # 表示 df_test
ダウンロードした時刻:2020/00/00 00:00:00 | |||||||
---|---|---|---|---|---|---|---|
NaN | 諏訪 | 諏訪 | 諏訪 | 諏訪 | 諏訪 | 諏訪 | 諏訪 |
年月 | 平均気温(℃) | 平均気温(℃) | 平均気温(℃) | 降水量の合計(mm) | 降水量の合計(mm) | 降水量の合計(mm) | 降水量の合計(mm) |
NaN | NaN | 品質情報 | 均質番号 | NaN | 現象なし情報 | 品質情報 | 均質番号 |
1970/1 | -3.3 | 8 | 1 | 49.5 | 0 | 8 | 1 |
1970/2 | -0.9 | 8 | 1 | 77.5 | 0 | 8 | 1 |
... | ... | ... | ... | ... | ... | ... | ... |
2019/9 | 21.3 | 8 | 1 | 29.5 | 0 | 8 | 1 |
2019/10 | 15.2 | 8 | 1 | 266.5 | 0 | 8 | 1 |
2019/11 | 7.7 | 8 | 1 | 23.5 | 0 | 8 | 1 |
2019/12 | 3.3 | 8 | 1 | 51.0 | 0 | 8 | 1 |
2020/1 | 1.4 | 4 | 1 | 22.0 | 0 | 4 | 1 |
604 rows × 1 columns
ん~一撃必殺。断然楽チンだ。圧倒的勝利。
「データファイルの読み込みにはPandasを使うべし!」
です。
Pandasでデータの形状を確認、整理
# 形状確認
df_test.shape
(604, 1)
何列かあるように見えて1列しかない。
# 統計的データ
df_test.describe()
ダウンロードした時刻:2020/00/00 00:00:00 | |
---|---|
count | 604 |
unique | 4 |
top | 1 |
freq | 601 |
先頭行がヘッダーとして読み込まれているので、本来の列が活かされない。
MultiIndex・ DataColumnsが1の時の対策
この状況への対応はどうすればよかったんでしたっけ。
- header=None : 新たなヘッダー行を自動追記の設定。
これはParserErrorが出る - header=0 : デフォルト設定と同じ。
状況変わらず。 - header=2 : 3行目(インデックス1番)をヘッダーにする。
今回はこれで上手く読めた(結果的に)
import pandas as pd # 気象庁よりDLした長野県諏訪のデータ data_file = "data.csv" # Pandasでcsvファイルの読込みD.F化 df = pd.read_csv( data_file, encoding = "shiftjis", header=2, #3行目をヘッダーにする index_col = 0 #1列目をインデックスにする ) # 表示(省略) #df
# 形状確認
df.shape
(602, 7)
# 情報確認
df.info()
<class 'pandas.core.frame.DataFrame'>
Index: 602 entries, nan to 2020/1
Data columns (total 7 columns):
平均気温(℃) 601 non-null float64
平均気温(℃).1 602 non-null object
平均気温(℃).2 602 non-null object
降水量の合計(mm) 601 non-null float64
降水量の合計(mm).1 602 non-null object
降水量の合計(mm).2 602 non-null object
降水量の合計(mm).3 602 non-null object
dtypes: float64(2), object(5)
memory usage: 37.6+ KB
使う行と列だけ抽出して上書きで置換してしまいます。
df = df.iloc[1:,[0,3]]
df
平均気温(℃) | 降水量の合計(mm) | |
---|---|---|
年月 | ||
1970/1 | -3.3 | 49.5 |
1970/2 | -0.9 | 77.5 |
1970/3 | -0.8 | 36.5 |
1970/4 | 9.1 | 71.0 |
1970/5 | 15.8 | 62.5 |
... | ... | ... |
2019/9 | 21.3 | 29.5 |
2019/10 | 15.2 | 266.5 |
2019/11 | 7.7 | 23.5 |
2019/12 | 3.3 | 51.0 |
2020/1 | 1.4 | 22.0 |
601 rows × 2 columns
# 欠損値の数量確認
df.isnull().sum()
平均気温(℃) 0
降水量の合計(mm) 0
dtype: int64
準備は整いました。
xとyの用意
グラフ作成に関して、x軸を日付、y軸に気温と降水量を個別に表したい。
Y軸にとる変数はy_kion
とy_kousui
って名前にしとこ。
X軸に日付のstr型で並べるとガチャガチャして見苦しいので、経過日数の数値に置き換えることにする。何月かわからなくなるけど。何か見た目に上下揺れてるのがわかりゃ良い。
配列作成や回帰分析をするので、DataFrameをNumpy配列に変換し作業する。
X軸の変数x作成
DataFrameをNumpy配列に変換する。
import numpy as np # dataframeをndarray変換 data_np = np.array(df) # またはdf.valuesでndarrayに #配列の形状確認 print(data_np.shape) # 要素の型 print(data_np.dtype)
(601, 2)
float64
行が600ヶ月、列が気温と降水量の2つを示す。
# 経過日数の配列xに変換 x = np.arange(data_np.shape[0]) + 1 # 配列表示 #x
Y軸の2変数作成
気温の配列 : y_kion
降水量の配列 : y_kousui
# 気温の配列 y_kion = np.array(data_np[: , 0]) # 気温配列の表示 #y_kion
# 降水量の配列 y_kousui = np.array(data_np[: , 1]) # 降水量配列の表示 #y_kousui
xとyの配列ができた。
グラフ描画
どんな感じかグラフで表示させてみる。
気温
# グラフ描画ライブラリのインポート import matplotlib.pyplot as plt # 気温のグラフ plt.plot(x, y_kion, color="red") # --- グラフ装飾(練習のため無駄に書く) --- # x軸ラベル名 plt.xlabel("Days") # y軸ラベル名 plt.ylabel("Temp") # タイトルラベル名 plt.title("Temperature\n1970/01 - 2020/01") # グリッド線 plt.grid() # ----------------------- # グラフ画像の保存 plt.savefig("temperature.png") # グラフの表示 plt.show()
おー。なるほどww
降水量
# グラフ描画ライブラリのインポート import matplotlib.pyplot as plt # 降水量のグラフ plt.plot(x, y_kousui, color="lightblue") # --- グラフ装飾(練習のため無駄に書く) --- # x軸ラベル名 plt.xlabel("Days") # y軸ラベル名 plt.ylabel("Rainfall") # タイトルラベル名 plt.title("Precipitation 1970/01 - 2020/01") # グリッド線 plt.grid() # ----------------------- # グラフ画像の保存 plt.savefig("precipitation.png") # グラフの表示 plt.show()
ははぁーん、なるほど。
Numpyで回帰分析
ここからnp.polyfit()
とnp.poly1d()
を使って回帰分析をしていきます。
# --- 方程式の算出 --- # 気温フィッティング 1 p_kion_1 = np.polyfit(x, y_kion ,1) p_kion_1
array([2.95345520e-03, 1.01179983e+01])
傾きと切片が得られた。
# 気温 f_kion_p1 = np.poly1d(p_kion_1)(x) #f_kion_p1
np.poly1d()
に、傾きと切片を持つnp.polyfit()
を渡すと回帰直線の方程式が得られるっぽい。
気温の回帰
# --- 描画 --- import matplotlib.pyplot as plt # リアルなデータのプロット plt.plot(x, y_kion ,label="Temp", color="red") # 回帰線 plt.plot(x, f_kion_p1, color="green", label="reg") # --- グラフ装飾(練習のため無駄に書く) --- # x軸ラベル名 plt.xlabel("Days") # y軸ラベル名 plt.ylabel("Temp") # タイトルラベル名 plt.title("Temperature\n1970/01 - 2020/01") # グリッド線 plt.grid() # 凡例 plt.legend() # ----------------------- # グラフ画像の保存 plt.savefig("temperature_2.png") # グラフの表示 plt.show()
1970年から2020年の50年間で平均気温が2~3度ほど上がってると見て良いのかこれは。
降水量の回帰
# グラフ描画ライブラリのインポート import matplotlib.pyplot as plt # --- 方程式の算出 --- # 降水量フィッティング 1 p_kousui_1 = np.polyfit(x, y_kousui ,1) #p_kousui_1
# 降水量 f_kousui_p1 = np.poly1d(p_kousui_1)(x) #f_kousui_p1
# --- 描画 --- # リアルなデータのプロット plt.plot(x, y_kousui ,label="Rainfall", color="lightblue") # 回帰線 plt.plot(x, f_kousui_p1, color="green", label="reg") # --- グラフ装飾(練習のため無駄に書く) --- # x軸ラベル名 plt.xlabel("Days") # y軸ラベル名 plt.ylabel("Rainfall") # タイトルラベル名 plt.title("Precipitation\n1970/01 - 2020/01") # グリッド線 plt.grid() # 凡例 plt.legend() # ----------------------- # グラフ画像の保存 plt.savefig("precipitation_2.png") # グラフの表示 plt.show()
多項式の次数を上げてみる
50次元まで。
# グラフ描画ライブラリのインポート import matplotlib.pyplot as plt # --- 方程式の算出 --- # 気温フィッティング 10 p50 = np.polyfit(x, y_kion ,50) #p50 # 気温 f_kion_p50 = np.poly1d(p50)(x) #f_kion_p100 # --- 描画 --- # リアルなデータのプロット plt.plot(x, y_kion ,label="Temp", color="red") # 回帰線 plt.plot(x, f_kion_p50, color="green", label="reg") # --- グラフ装飾(練習のため無駄に書く) --- # x軸ラベル名 plt.xlabel("Days") # y軸ラベル名 plt.ylabel("Temp") # タイトルラベル名 plt.title("Temperature\n1970/01 - 2020/01") # グリッド線 plt.grid() # 凡例 plt.legend() # ----------------------- # グラフ画像の保存 plt.savefig("temperature_3.png") # グラフの表示 plt.show()
/data/data/com.termux/files/home/py3.8/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3319: RankWarning: Polyfit may be poorly conditioned
exec(code_obj, self.user_global_ns, self.user_ns)
意味無さげなふにゃふにゃ線。