気象データアナリスト講座PBL

コインランドリーの利用者数と気象との関係を分析する

本PBLではコインランドリーの稼働履歴データと気象データを使い、気象条件がコインランドリー利用者数に与える影響を分析します。
またその結果をもとに、気象条件によって起こる稼働のピークを平滑化する施策を立案します。

分析には、ウォッシュプラス富士見1丁目店(千葉県浦安市)のデータを利用します。

分析の流れ

  1. 稼働履歴データを集計して稼働回数を求め、稼働回数を利用者数とみなして考察します。
  2. 気象データを稼働回数データに結合し、気象条件による稼働回数の変化を分析します。
  3. 稼働回数を気象データ等でモデリングし、気象の影響を多角的に分析するとともに、明日の稼働回数の予測を試みます。

コインランドリーの利用者数と気象との関係の仮説


目次

1. 稼働履歴データの集計・可視化

1-1. 稼働回数・売り上げの集計

まず必要なライブラリをまとめてインポートします。
ここで必要に応じて、グラフ描画ライブラリmatplotlibの日本語表示設定も行います。

washplusフォルダにある、稼働履歴データを読み込みます。
稼働履歴データのファイルは通常のCSV形式なので、Pandasのread_csv関数で簡単に読み込めますが、ファイル数が複数あります。

ウォーミングアップとしてまずファイルを1つ、どれでも良いので読み込んでみましょう。
ポイントとして文字コードがshift-jisのため、read_csv関数のオプションとして指定する必要があります。
読み込んだら最初の5行を表示させてみましょう。

次に全ファイルをまとめて読み込みます。
ここではglob関数を使って全ファイル名を取得したのち、forループで順番に読み込んで、concat関数で縦につなげる処理を行います。
こうすることで、全データを1つの変数に格納します。

読み込んだファイル名を見ると、日付順になっていないかもしれません。
その場合、変数laundry_orgの中身も日付順になっていませんので、日付の古い順に並び替える処理をします。

ここでは、Pandasのto_datetime関数で開始日時をTimestamp型に変換します。
Timestamp型にすると、年・月・日・時・曜日などの情報に分解でき、また日付の足し算(e.g. 1時間後を計算)などができて便利です。
Timestamp型に変換したら、sort_values関数で開始日時の古い順(降順)にソートします。

2018年2月〜2020年9月の稼働履歴データを読み込みました。

さっそく幾つかの切り口で集計をしてみましょう。
カラム名を見ると、機器種別・号機・コース番号などで集計できそうです。また稼働回数だけでなく、売上の計算もできそうです。

機器種別・号機・コース番号について
|号機|機械種別| |--|--| |1|乾燥機(二段)| |2|乾燥機(二段)| |3|乾燥機(二段)| |4|乾燥機(二段)| |5|洗濯乾燥機(中型)| |6|洗濯乾燥機(中型)| |7|洗濯乾燥機(中型)| |8|洗濯乾燥機(大型)| |9|洗濯乾燥機(大型)| |10|洗濯乾燥機(特大)| |コース番号|コース名| |--|--| |1|洗濯乾燥(標準)| |2|洗濯乾燥(少なめ)| |3|洗濯| |4|乾燥|

各カラムについては付属資料「washplusデータ基本情報.xlsx」も参照して下さい。

ここでは気象との関連が大きそうなコース番号を中心に見ていきます。

まず稼働回数から集計します。
value_counts関数を使うと、指定したカラム(e.g. 機器種別)の構成要素(e.g. 洗濯乾燥機, 乾燥機)について出現回数をカウントできます。

乾燥コースの利用頻度が最も多く、次いで洗濯乾燥(標準)コースが利用されています。


次に支払料金カラムを使って、売上を計算します。

groupby関数を使うと、指定したカラム(e.g. 機器種別)の構成要素(e.g. 洗濯乾燥機, 乾燥機)でグルーピングし、集計したいカラム(e.g. 支払料金)についてグルーピングした構成要素ごとの合計(sum関数)を計算することができます。

稼働回数とは逆に、洗濯乾燥(標準)コースの方が売上は大きくなっています。
稼働回数は乾燥コースより少なくても、単価が大きい洗濯乾燥(標準)コースの方が売上は大きくなるようです。

目次へ戻る

1-2. 月ごと・コースごとの稼働回数の集計

ざっくりと集計してみましたので、今度は日付情報を使って月ごとや週ごとに利用者数(稼働回数)の特徴を見ていきます。

事前準備として、日付情報を年・月・日・時・曜日に分けて格納します。
上で開始日時をTimestamp型に変換してあるので、年・月・日・時・曜日の情報は以下のようにして簡単に取り出せます。

また月ごとに集計する都合上、ちょうど2年分のデータに絞ります。今回は2018年9月〜2020年8月を使います。

今後、何回か積み上げ棒グラフで可視化するため、まず関数を定義しておきます。

次に、月ごと・コースごとにクロス集計し、各コースの稼働回数が月(季節)によりどう変化するのか見てみます。

crosstab関数を使うとExcelのピボットテーブルのような集計(出現回数のカウント)ができます。

以下のセルでは、クロス集計結果を上で定義したbar_plot関数に与えてグラフを作成し、クロス集計結果そのものも出力させています。

乾燥コースは7月が突出して多いことがわかります。ついで12月が多くなっています。
洗濯乾燥(標準)コースは乾燥コースに比べると月ごとの差は小さいものの、12月・7月の順に多くなっています。

これらの月で稼働回数が多いのは何故でしょうか?
7月は梅雨、12月は年末であることが関係していそうですね。

目次へ戻る

1-3. 曜日ごと・コースごとの稼働回数の集計

続いて曜日ごと・コースごとにクロス集計し、曜日による利用者数の変化を見てみましょう。

各コースとも土日、特に日曜日の稼働回数(利用者数)が多いことがわかります。

目次へ戻る

1-4. 時間ごと・コースごとの稼働回数の集計

今度は時間ごと・コースごとにクロス集計し、時間帯による利用者数の変化を見てみましょう。

全体として午前中に増加する稼働回数は10時ごろにピークに達し、夕方ごろまで多い状態が続いたのち、19時ごろから減少する傾向があります。
コースごとに見ても、洗濯乾燥(標準)コースと、乾燥コースは全体と同じような傾向があります。
一方で洗濯乾燥(少なめ)コースに関しては、午前10時ごろよりも夕方の方が稼働回数が多くなっています。


ここまでの分析結果を簡単にまとめます。


目次へ戻る

2. 気象データを結合して気象条件で集計・可視化

2-1. 天気とコースでクロス集計

ここからは稼働履歴データに気象データを結合して、気象状況の切り口でデータ分析を行っていきます。

分析は基本的に日単位に集計した稼働回数と、日単位に集計した気象データで行います。

まず稼働回数を日単位に集計したデータフレームを作るため、crosstab関数を使ってdateとコース番号でクロス集計を行います。

次に気象データを読み込みます。気象データはweatherフォルダの中に格納されています。
今回はアメダス江戸川臨海の観測データを、気象庁ホームページより事前にダウンロードしておいたものを利用します。

ダウンロードしたファイルはCSV形式ですが、ヘッダが最初の行(0行目)にないなど、特殊なフォーマットとなっています。
そのため、うまく読み込むためには工夫が必要です。

ここでは以下のように、アメダスCSVデータを読み込むための関数を定義して使用します。

ダウンロードしたアメダスCSVにはカラム名として使える適切な行がないため、あらかじめ用意したカラム名リストを渡して読み込みます。
それが引数all_paramsです。

事前に用意したアメダスCSVデータは、品質情報と均質番号を含んだ状態で入手しました。
今回は品質情報が8または5の場合は良しとしますが、5未満の場合は欠損値として扱うことにします。
また均質番号は期間通して1であることを確認済みとし、今回は使用しません。
品質番号はカラム名を気象要素名_q、均質番号は気象要素名_iとすることで識別できるようにカラム名を付けます。

読み込んだデータのうち、均質情報などは読み込み時にしか使いませんので、関数からreturnする際に落とすことにします。
必要なカラム名のみのリストを渡してその処理をします。それが引数paramsです。

それではアメダス江戸川臨海の日別観測値を読み込みましょう。

アメダス江戸川臨海は4要素(気温・風向風速・降水量・日照時間)の観測はありますが、相対湿度の観測がありません。
しかしながら、コインランドリーと相対湿度には何らかの相関があることも考えられ、できれば相対湿度も使いたいです。
そこで東京管区気象台の観測データから相対湿度だけダウンロードしました。これを利用することにします。

アメダス江戸川臨海と東京管区気象台のデータを結合します。
merge関数を使い、年月日カラムを軸にして、東京管区気象台のデータをアメダス江戸川臨海に結合します。

これで日別の観測値を読み込むことができました。しかし、これだけでは晴れ・曇り・雨の天気カテゴリーが不明です。
東京管区気象台の天気の観測を使う方法もありますが、今回はせっかくアメダス江戸川臨海が近くにあるので、日照時間や降水量の観測から天気カテゴリーを推定する方法をとります。
(なお天気カテゴリーは、推計気象分布から入手する方法もあります。)

以下のようなロジックで天気を推定します。

  1. 降水量があれば雨とする(雪は考えない)
  2. 降水量がゼロの場合、日照時間で晴れか曇りを判定する
    • あらかじめ月毎の日照時間の最大値を求めておく(毎月一日くらいは日中ずっと日照がある日があると仮定)
    • その最大値と、日毎の日照時間の比を計算し、0.2未満であれば曇り、それ以外を晴れとする

(注意)
この方法で天気を推定した場合、例えば「早朝まで雨が降っていたが日中は晴れた」ケースであっても日天気が「雨」と判定されてしまう

ではまず観測日(年月日)をTimestamp型にして月を抽出し、月でグルーピングし、月ごとの日照時間の最大値を求めます。

概ね違和感のない値が求まりましたので、monthカラムを軸として、変数max_sun_hourを変数obs_dailyにマージします。
そうすると、変数obs_dailyに日照時間maxカラムが結合されます。
次に日照時間と日照時間maxの比を計算し、日照時間rateというカラムを作ります。

日照時間rateと降水量から天気カテゴリーを作ります。

これで日単位のコインランドリー稼働回数データと気象データが準備できましたので、dateカラムを軸に両者を結合します。

これで天気とコースごとの分析を行う準備ができました。

ではさっそく天気とコースごとにクロス集計し、稼働回数の平均値・標準偏差・中央値を求めてみましょう。

天気が悪くなるほど、乾燥コースの稼働回数が増えています。
その他のコースはあまり天気で差がありませんが、あえていうと洗濯乾燥(標準)コースは雨の日に少しだけ少ないようにも見えます。

また稼働回数が多いほどばらつきも大きくなっているようです。この様子を箱ひげ図で可視化してみます。

目次へ戻る

2-2. 日照時間とコースごとの分析

天気カテゴリーによる分析では、雨天ほど利用者数(稼働回数)が多い傾向がありました。
そうすると日照時間と利用者数には逆相関があるかもしれません(日照時間が増えるほど稼働回数は減る)。

そこで、今度は連続値の変数である日照時間と利用者数との相関係数を計算してみましょう。Pandasのcorr関数を使います。
なお、これより先は稼働回数が多い乾燥コースと、洗濯乾燥(標準)コースに絞って見ていきます。

日照時間とコース1はほとんど相関がありませんが、コース4は負の相関係数が0.5以上あり、それなりに相関がありそうです。

次に、日照時間と稼働回数はどちらも連続値のため、散布図を描いてばらつき具合を可視化してみましょう。

→ 天気カテゴリーで得られた結果と同様の傾向が見られました。

目次へ戻る

2-3. 連続降水日数とコースでクロス集計

今日の天気だけでなく、昨日・一昨日の天気も利用者数に影響を与えそうです。
例えば雨の日が数日連続して続くと、自宅で干せないため、コインランドリーで洗濯や乾燥をする人が増える、という仮説が立てられます。

そこで、当日を含めて雨が何日続いているか、すなわち連続降水日数をもとに分析を行います。

連続降水日数の求め方は、天気カラムの値をshift関数で1日ずつずらしたカラムを作り、何日連続でRainyの天気が続いているかを逐次計算していく方法をとっています。

どうも梅雨の時期に10日以上連続して雨天が続いたケースがあるようです。
今回はサンプル数のバランスもふまえて、1日連続と2日連続を1つにまとめます。値としては2日連続の場合に1を代入します。
また3日以上連続も1つにまとめます。値としては4日以上連続の場合に3を代入します。
これで連続降水日数を、0日連続|1〜2日連続|3日以上連続、の3つのカテゴリーに分けたことになります。

連続降水日数が求まりましたので、改めてコインランドリー稼働回数データと気象データを結合します。

では連続降水日数とコースごとにクロス集計しましょう。
統計量としては中央値を使い、またコースも洗濯乾燥コース(標準)と乾燥コースに絞って見ていきます。

目次へ戻る

2-4. 降水量とコースごとの分析

降水量が増えると利用者数にどのような影響があるでしょうか?

乾燥コースの利用者数が増えそうですが、逆にしっかり雨が降っていると、そもそも外出を避ける傾向もあるかもしれません。
そこで、まず降水量とコースごとの利用者数との相関係数を計算してみます。

降水量と乾燥コースの相関係数はプラスの値になりましたが、相関係数自体は大きくありません。
次に、降水量と稼働回数はどちらも連続値のため、散布図を描いてばらつき具合を可視化してみましょう。

目次へ戻る

2-5. 相対湿度とコースごとの分析

湿度が上がると利用者数にどのような影響があるでしょうか?

湿度が高いと洗濯物が乾きにくいため乾燥コースの利用者が増えそうです。
(ただし、単純に雨天の時は湿度が高いので、雨の影響を間接的に見ていることになるかもしれません。)

まず相対湿度とコースごとの利用者数との相関係数を計算してみます。

湿度と洗濯乾燥コースには相関がなさそうですが、乾燥コースは相関係数0.5以上あり、それなりに相関がありそうです。
また平均湿度より最小相対湿度の方が少しだけ相関係数が大きくなっています。

そこで最小相対湿度について、散布図を描いてばらつき具合を可視化してみます。

→ ただし、雨天の時に湿度が上がる様子が現れているだけかもしれないので、解釈には注意が必要です。

目次へ戻る

2-6. 水蒸気量とコースごとの分析

先ほどは空気の水分量を相対湿度で扱いましたが、空調の分野では絶対湿度、すなわち水蒸気量も扱うようです。
確かに同じ相対湿度でも、気温によって空気中の水蒸気量は異なります。
洗濯物を干す場合でも、空気があとどれくらいの水蒸気量を含むことができるか、という指標が洗濯物の乾き具合に影響しそうです。
つまり、空気が追加で含みうる水蒸気量が小さいと、洗濯物が乾きにくくて乾燥コースの利用者が増える、との仮説を立てられます。

そこで、平均気温と平均湿度から飽和水蒸気量と水蒸気量を計算し、両者の差分を取ることで「追加で含みうる水蒸気量」を求めてみます。
計算にはwxparamsを利用します。
また「追加で含みうる水蒸気量」では長いので「可蒸発量」という造語でカラム名をつけることとします。

想定した通り、乾燥コースの稼働回数と可蒸発量には、それなりの相関がありそうです(負の相関)。
また洗濯乾燥(標準)コースと稼働回数には相関があまりなさそうです。

では散布図も描いて、ばらつき具合を可視化してみましょう。

目次へ戻る

2-7. 午前中の天気とコースでクロス集計

先ほどは単純に日積算の降水量と日照時間で天気を推定しました。
しかしコインランドリー利用者、すなわち洗濯する人の行動や気持ちを想像すると、午前中の気象条件の方が影響が大きいかもしれません。

そこで今度は時別の気象データから午前中の天気を推測してみます。

まず時別のアメダス観測値を読み込みます。
時別のアメダス観測値は複数のファイルに分かれているため、まずglob関数でファイル名のリストを取得します。
そして日別アメダス観測値を読み込む際に作ったRead_AMeDAS_CSV関数で読み込みを行います。

これで時別のアメダス観測値を読み込むことができました。

次は午前中の天気の推定です。
日天気の推定と同様に、午前中の日照時間rateと降水量を使って天気を推定します。

午前中の天気が求まりましたので、変数obs_dailyに結合します。
また稼働回数データにも結合します。

では、これまでと同様に午前中の天気とコースごとにクロス集計しましょう。
統計量としては中央値を使い、またコースも洗濯乾燥コース(標準)と乾燥コースに絞って見ていきます。

目次へ戻る

2-8. ゲリラ豪雨(夕立)の有無とコースでクロス集計

朝は晴れていたので洗濯物を干して出かけたが、急なゲリラ豪雨(夕立)で洗濯物が全滅…
なんてことがあると、乾燥のためコインランドリーを利用する、という人が増えそうです。

そこで今度は時別の気象データから夕立の有無を判別し、利用者数の違いを分析します。

夕立を次のように定義します。

夕立の定義を満たす日は、12日しかないようです。
ちょっと少ない気もしますが、今までと同様に夕立の有無で集計したり、箱ひげ図を描いてみます。

以下の条件の日を抽出して比較を行います。

目次へ戻る


ここまでの分析結果を簡単にまとめます。


目次へ戻る

3. 各コースの1日あたり稼働回数を気象データ等でモデリング

3-1. 仮説検定と統計モデリング(線形回帰)

これまでは主にクロス集計と可視化を用いて、コインランドリーの利用者数(稼働回数)と気象状況との関係性を分析してきました。
この節ではより定量的に稼働回数と気象データとの関係性を見ていくため、統計ソフト『R』を用いて作業します。

これまでに整形したデータフレームを、Rで使うためにCSVファイルで出力します。


→ Rで仮説検定と統計モデリング(線形回帰)へ


目次へ戻る

3-2. 機械学習で予測(線形回帰・ランダムフォレスト)

3-2-1. 線形回帰

今度は機械学習で予測にチャレンジします。
Rでも使った線形回帰と、機械学習の代表的なアルゴリズムであるランダムフォレストを使ってみます。

予測ターゲットはより気象との関係性が高かったコース4の稼働回数とします。

特徴量は、Rで行った変数選択と同じ特徴量を採用します。

Pythonの場合、カテゴリー変数はダミー変数化するなどの前処理が必要になります。
ここでは連続降水日数と月もカテゴリー変数として扱うため、型を文字列に変換します。
そしてカテゴリー変数をダミー変数化する処理をします。

ダミー変数化したら、多重共線性を回避するため、各カテゴリー変数とも1種類のカテゴリーを削除する必要があります。
ここでは、Rの結果に倣ってカテゴリーの1種類を削除します。

機械学習では、データを学習用・検証用・テスト用に分けて予測モデル開発を行います。
ここでは最初の1年分を学習用、残りの1年分を検証用とします。
(またテスト用データについては後述します。)

ちょうど同じデータ数に分割されました。

次に予測ターゲット(y)と、特徴量(X)に分割します。

データの準備が整いましたので、線形回帰から試してみましょう。

余談ですが、機械学習はただ使うだけなら数行のコードでできてしまいます。
機械学習を使う前までの、データの処理がいかに大変かがわかるかと思います。

目次へ戻る


3-2-2. ランダムフォレスト

それでは次にランダムフォレストを試してみましょう。
ランダムフォレストの場合、多重共線性を気にする必要はないので、全気象要素を特徴量にしても良いのですが、線形回帰との比較のため、また簡単のため、線形回帰と同じ特徴量で学習を行うこととします。

一般にはランダムフォレストの方が高精度を期待できますが、どうなるでしょうか?

ランダムフォレストはハイパーパラメータの調整で予測精度が変わりますが、今回は木の数(n_estimators)を500本にしてみました。
また学習の再現性を確保するため、randam_stateを固定しています。

結果を見ると、MAE・RMSEともに線形回帰の方が精度が良い結果となっています。
どのような予測が出ていたのか、散布図で可視化してみます。

線形回帰に比べると、ランダムフォレストの方が大外しをするケースが多いため、RMSEが悪くなったと考えられます。
またランダムフォレストはアルゴリズムの特性上、実績の下限値付近を予測しにくい傾向が現れています。

今回、ランダムフォレストに入力した特徴量は線形回帰と同じものにしました。
しかし、ランダムフォレストの場合は多重共線性は気にしなくて良いので、入力できる全ての特徴量を利用しても良いのかもしれません。
ぜひ試してみて下さい。

では、各モデルがどのような学習をしたか見てみましょう。
まず線形回帰の回帰係数をデータフレームで書き出します。

Rの線形回帰と全く同じにはなりませんが、同じような傾向(正負の符号)で回帰係数が求まっていると思います。

次にランダムフォレストで重要度を見てみましょう。

重要度の上位には、最小相対湿度や日照時間など天気を示唆する特徴量や、休日・前日の実績のような非気象の特徴量が入っています。

目次へ戻る


3-2-3. GPVを使った予測

ここまでで、機械学習のモデル開発ステップのうち、学習と検証を行いました。
残るテストは、数値予報モデルMSMからある1日の予測値を切り出して学習済みモデルに投入し、予測を出してみます。

MSM-GPVの処理には、気象ビジネス推進コンソーシアムのwxbcgribを利用します。
また一部の気象要素の処理にwxparamsを利用します。

今回は2020年9月6日のコース4の稼働回数を、MSM 2020年9月5日12zベースを使って予測します。
data_20200906.png

MSM(GRIB2形式)は気象フォルダに入っています。
wxbcgribを使い、ウォッシュプラス富士見1丁目店の場所の気象予測データを取得します。

予測対象である9月6日の予測データに絞ります。

これで9月6日対象の1時間ごとの予測データを取得できましたが、まだ予測はできません。
学習済みモデルを作った時と同じ特徴量(1日ごとの気象要素など)に変換する必要があります。

今回は1日分のテストデータを作るだけなので、少々強引に1つずつ特徴量を計算していきます。
まずテストデータを入れる「箱」を用意します。

最高気温・最大風速・最小相対湿度は、それぞれ1時間ごとの値の最大値・最小値を計算して求めることができます。

日照時間をどのように求めるかは工夫が必要です。
日射量を使う手もあるかもしれませんが、前1時間の平均値から求めるのは困難です。
ここでは全雲量を使い (100 - 全雲量) / 100 を日照時間の推定値として使用することとします。

午前中の天気と降水量です。
まず数値予報は非常に小さい値の降水量を計算する場合があるので、0.1mm未満の微弱な降水量は切り捨ててゼロとします。
その後、午前7時〜午前12時の予測データを用いて、午前中の天気と降水量を計算します。

連続降水日数は、前日9月5日までの連続降水日数を使い、9月6日の予想天気が雨ならば+1することで求まります。
上で示した表を見れば、9月5日までの連続降水日数がわかります。

残りは前日9月5日の稼働回数の実績を代入し、9月と日曜日のカラムを1にします。

やっとテストデータが完成しました!
では2つの学習済みモデルに投入して、予測を出力しましょう。

線形回帰・ランダムフォレストともに、良い予測ができたと言えるのではないでしょうか。


本PBLでは気象観測データを使って学習モデルを作り、気象予測データを入力してターゲットの予測を得ました。
このようなアプローチをPPM (Perfect Prognosis Method)と呼びます。

これに対し学習モデルを作る段階から過去の気象予測データを使い、気象予測データを入力してターゲットの予測を得る方法もあります。
このようなアプローチをMOS (Model Output Statistics)と呼びます。

気象データは同じ気象要素でも、観測データと予測データという性質の異なる2種類のデータが存在します。
観測データと予測データ、PPMとMOSの違いを理解して使い分ける(使いこなす)ことが重要です。


これまで、コインランドリーの利用者数と気象状況との関係性を分析し、様々な知見を得るとともに、予測まで行うことができました。
これらの結果をもとにして、三密回避&機会ロス削減のためのピークシフト施策を考えてみましょう。

目次へ戻る