専門実践教育訓練講座の最終成果物として記事作成。
Kaggleの上位3%以下という結果。
かなり冗長となるが講座の合否に使用される為、分析順序と考察を交えて記載。
目次
- はじめに～記事概略まとめ
- 本記事について
- 分析環境
- テーマについて
- データ分析1：データ確認
- 最初の機械学習
- 説明変数の編集
- 目的変数の扱い
- L2正則化（Ridge回帰）
- 機械学習2
- LightGBMによる予測
- スコアを上げるための施策
- まとめ
はじめに～記事概略まとめ
- Kaggleコンペ、Bike Sharing Demandのデータ分析を実施
- スコアは上位3%に収まる
- Aidemy様のデータ分析講座を受けた最終成果物として本記事作成
- 本記事では考察を重視しており、かなり詳細まで書いている（＝冗長）
具体的な実施項目。
- 内容：ワシントンD.C.に於けるレンタル自転車の利用量（1時間当たり）
- 目的：利用者数を予測。天候、日時、休祝日情報より
- 分析
- 各データについて確認
- 事実上の外れ値あり → 実際の天気予報を入手し補完
- 「体感温度」の外れ値はNWS（アメリカ国立気象局）紹介の理論式を参照し計算で補完
- データに施されているノーマライズを解読、再処理
- 線形理論での学習
- 線形モデルでは不適と判断
- LightGBM（決定木ベース）でモデル作成
- optunaによるハイパーパラメータの最適化
- 特徴量の再考
- 「ユーザ属性」情報がない → 模した特徴量を作成
- 時系列上で利用者が急峻に下がる点があり → この特徴量を作成
- 逆に、ピークとなる箇所の特徴を作成
- クリスマス以降を休日扱いに変更
- 目的変数の生起確率分布を確認、log(x+1)で変換することに
- LightGBMで再学習、評価結果は上位2.9%
- まとめ
- 各データについて確認
以下はこれらの詳細となる。
本記事について
作成の経緯
Aidemy様のデータ分析講座の最終課題にて作成。
以下を満足すること、という指定あり。
- 再現性のある事
- 機械学習の手法（各処理の手順等）は講座で学習した内容をカバーしているか
- 考察・反省はあるか
著者について
データ分析歴
- 機械学習・ビッグデータという「データ分析」は未経験
- 物理実験およびそのデータ分析に関する経験は豊富
プログラミング経験
- 業務でWinWrap Basic（商材の実験機器内蔵Editor）での機器制御、信号処理、科学演算、解析等
- 学生時代にROOT（CERNによるデータ開発環境。C++ベース）で実験データ分析
- 業務自動化やweb予約自動化等でPython 3 デスクトップアプリ作成
- 趣味でArduino（C++）を少々
プログラミングを専門とした業務は行っていないのでアマチュアレベル
分析環境
| 使用環境 | |
|---|---|
| Kaggle Notebook | 2024年5月現在の仕様 |
| 言語 | Python3 |
| Google Chrome | ウェブブラウザとして |
| Windows 11 | 使用OS |
テーマについて
Bike Sharing Demand 概略
Kaggleのコンペティション、Bike Sharing Demand。
ワシントンD.C.におけるレンタル自転車の利用者数の予測問題。
1時間ごとに利用者数と気象データ等の周辺データを記録したもの。
2011年1月1日～2012年12月31日までのデータ。
周辺データは、「日付」「天候」「風速」「湿度」「祝日か平日か」等。詳細は後述。
1時間当たりの利用者数が予測対象（目的変数）。
訓練用データは毎月1日～19日までのデータ、テスト用データは毎月20日以降のデータが与えられている。
予測対象は「各月の20日以降の利用者数」となる。
適当な図だが、予測対象を視覚的に説明。
テーマ選定理由
以下の点を考慮。
- マーケティングらしい題材
- ビッグデータらしい収集データ
- 一般に理解・イメージしやすい
- 説明変数と目的変数が少なく、分析・考察・説明しやすい
- 結果の評価をしやすい（Kaggleのコンペにて、順位が分かる）
分析用データ
Kaggle提供のデータ。
提供データはこちら（Kaggleのデータタブのリンク）。
データの出典はこちら。
見かけ上欠損値がなく、扱いやすいデータ。
Bike Sharingサービスの仕様上自動的にデータ取得される為か。
但し、ところどころ記録のない時刻（行ごと無し）があり、また事実上の欠損値となるデータがある。
また、データの説明書きと実際のデータとで形式が異なる点があることに気付く。
これは分析を進める中で触れる。
評価指標
評価はRMSLE値 (Root Mean Square Log Error の略)で、レンタル数の予測値と実際のレンタル実績数との差を評価。
RMSLE \ = \ \sqrt{\frac{1}{n}\sum^n_{i=1}\Bigl(log (p_i+1) - log(a_i +1)\Bigr)^2}
- n : データ数。1時間を1データとして扱う
- $p_i$： 予測値
- $a_i$： 実際の値
- $log$： 自然対数
目標値設定
上位5%以上を目標に定める（Kaggleの「Leaderboard」上の順位）。
RMSLE値およそ0.39が上位5%。これを目標にする。
「RMSLE値がxxxの場合、モデルの精度がこの程度だ」といった数値からの精度目標設定は難解なので、
目標はRMSLE値からではなくKaggleの順位にした。
凡そ偏差値66%程度を目標にした。
目標RMSLE値について掘り下げ
「RMS」はよく使用するが「RMSLE」は今回初めて触れるので少し詳しく考える。
結局、図的に理解すると、
- RMSE ： 切片で誤差の幅を見る
- RMSLE ： 傾きで誤差の角度を見る
といったイメージ。
冗長になるのでRMSLEの詳細は別の記事にする。
| RMSE | RMSLE |
|---|---|
 |
 |
データ分析1：データ確認
データcsvは、Kaggle Notebookの場合は以下よりダウンロードできる。
\kaggle\input\competision name*.csv
csv種類：train.csv、test.csv、sampleSubmission.csv
Kaggleでの一般的なものなので、詳細割愛。
各変数
| 変数名 | 内容 |
|---|---|
| datetime | 日時 |
| season | 1 =春、2 =夏、3 =秋、4 =冬 |
| holiday | 1=祝日、0=他（土日は0扱い） |
| workingday | 1=平日（土日祝を除く日）、0=土日祝 |
| weather | 1=晴れ、2＝曇り・霧、3＝軽い雨雪、4＝強い雨雷雪 |
| temp | 気温(℃)。アメリカだがセルシウス温度 |
| atemp | 体感気温(℃) |
| humidity | 湿度(%) |
| windspeed | 風速（単位不明） |
| casual | 利用者数、うち非会員 |
| registered | 利用者数、うち会員 |
| count | 利用者数 |
countが目的変数。
casual、registeredも説明変数の内訳という意味で目的変数に近い。
test.csvにはcount、casual、registeredは含まれない。
目的変数の配置は前述の通り。
情報要約
一部重複するがまとめておく。
- 欠損値なし（見かけ上。事実は一部欠損値は0または0相当に置換されている模様。後述）
- 説明変数は9種類
- 目的変数はcount、int型
- count、registered、casualはtrain.csvのみに含まれる
- データの日付は2011-1-1 00:00:00 ～ 2012-12-31 23:00:00の間
- train.csvは毎月1日～19日までのデータ
- test.csvは毎月20日以降のデータ
- datetime以外の説明変数は全て数値型（カテゴリカルなものはエンコード済み）
- datetime（日付に相当する説明変数）の型はobject型
- countは最小で1人、最大で977人と、かなり数値の開きがある
ちょっとした疑問：temp、atemp
気温、体感温度が0を下回らないようになっている。
マイナスになると計算値の符号が反転してしまうため、セルシウス温度は科学計算向きでないことはよく経験する。
今回はどのように補正されているのだろうか？
詳細調べると、Kaggleの以下のページに各変数の簡単な説明がある。
https://www.kaggle.com/competitions/bike-sharing-demand/data
また、Dataの出所を調べると以下ページが紹介されている。
https://archive.ics.uci.edu/dataset/275/bike+sharing+dataset
こちらを見ると、そもそもatemp、tempは「Normalizeした値」として説明されている。
ここでは値が0~1となるような式が書かれており、取得したデータと合致しない。
また、体感気温のatempについて。
気温、風速、湿度から計算されるものと思われる。
日本では一般にミスナールの計算式・改良版で体感温度を出すらしい。
しかし、temp、windspeed、humidityを使用して計算してもatempとは似つかない結果を得た。
この疑問から出発し、天気予報データを取得、温度情報のノーマライズ式を検証、デノーマライズして、体感温度を再計算、体感温度をノーマライズし、結果的にatempの欠損値を埋める、という処理を経るに至る。
詳細は別途、追って別記事にしたものを紹介する。
Kaggleやウェブ上で色々調べたが、この検証をした人はいないと思われる。
追って紹介するが、以下記事に別途まとめた。
https://qiita.com/SumMae/items/d0548ee2156965bdde12
train、testの結合：「Test_Flag」属性
Aidemy様で常套手段として。Test_Flagというカラムを作成。
データ分析時にtrainとtestに偏りがないか調べる。
説明変数がtrainとtestで異なると予測に使えない。
（特徴量エンジニアリングの際にtrain、test双方に処理をする必要がある）
これに便利な処理としてdf_allとしてtrainとtestを結合する。
結合後に分類が分かるように、「Test_Flag」とうカラム追加、=1ならtestデータ、=0ならtrainデータとする。
Test_Flag、という属性がしばしば出てくるが、こうして作られた属性である。
Test_Flagで分類した諸データを表示する。
weather=4はほとんど無い。予測するうえで統計的に数が少ないので、weather=3に含めてしまうことにする。
他、副次的に得られる情報。
windspeedに0が多く、欠損値の可能性がある。
humidity=0もあるが、地球上では有り得ない数値なので欠損値と断定。
年、月、日、時刻、曜日への分割
ほぼ前処理不要なデータが提供されているが、「datetime」だけ修正したい。
種々プロット・集計などして性質を可視化したいが、datetimeは年月日時刻の情報を1変数に集約した、全行に亘りユニークな値となっており時系列分析以外では扱い難い。
また、曜日情報も参考として欲しい。
曜日は、経験例として金曜昼はムスリムの礼拝で路上駐車が多くて通行に悩まされるという経験があり、「金曜昼だけ近くの例は異常にレンタル自転車で行く、という集団がいるのでは？」と期待。
結果、その特徴は顕在化しなかったのだが。
併せて、これらデータは最終的にcategory変数属性にしておく。
（ただし、ソートや特徴量作成時の範囲設定時には一時的にint型として扱うことがある。
category型の場合、「["hour"]<12:00」といった範囲条件抽出ができない為）
年月日等データ確認
2011年、2012年比を出したいことと、コードが正しいか確認目的もあり、総数での表示。
各月のcount数を表示。
各月の傾向はおおよそ同じ。
- 2011年count数合計：781,979
- 2012年count数合計：1,303,497
- 2012年/2011年比：約1.67
seasonでも同様。
season{1:1~3月、2:4~6月、3:7~9月、4:10~12月} という括りであること確認。
monthとseasonとは同じ意味で分割するビンの違い、というイメージ。
曜日、時刻、workingday、holiday関連。
特に平日と休日とでは分布が大きく異なる。
天候データ確認
weatherについて。
countの総数表示。
各天候がいくつ出現したか表示。weather=4は3回しか出現していない。
ratioは切り捨てられて0%扱いになるくらい出現頻度低い。
天候は1が圧倒的に多く、count値に於ける71%を、全行に占める割合（出現回数）で66%が天候=1。
天候=1,2,3,4で降順となっている。
4はほとんど見られない。
各天候に於けるcount値の平均値を時刻に対してプロットすると、分布の形はほぼ同じで量が小さくなる。
この分布の形状は先に示した、全日付に対して合計した「count vs hour」の形状にほぼ一致している。天候による変化は単純な全数の変化のみで、どこかの時刻等に対する特別な変化はしない模様。
casual、registeredについて確認
count = casual + registered に分かれる。
これらの確認。
特筆すべきは、ヒストグラム。
casualは横軸の絶対値を無視するとλ=1のポアソン分布のような形状で、レンタル数は正規分布のような変動ではなくゼロを起点として右裾長な分布を示す。
registeredも似ているが、100あたりにも緩いピークを持つような分布。
他、casualでは特に季節的に気候が良い（monthに於ける4月～9月）と利用者数が増えてcountに占める割合も上がる。
量的変数の可視化
seabornライブラリで一度に描いた。
- temp、atempは相関が0.99、ほぼ傾き=1の直線上に乗る
- casualはregisteredと比較してtemp、atempが低いと数が顕著に少なくなる
- humidity、windspeedは全体の数値から離れた位置に0となる値を示す。欠損（データ取得できなかったケース）を0で埋めている可能性が高い
- humidityが高い個所でcasualの分布が顕著に減る
- windspeedが高い個所でcasual、registered、countそれぞれの分布が減る
- casualとregister（またはcount）間の分布をみるとハート形となっており、いわば2種類の異なる系統の確率分布が存在しているような分布の仕方を示す
余談だが、atemp=30あたりでhumidity、casualの分布が割れており、またatemp=20あたりでwinspeedの分布が割れている。
atempとtempとは相関の強い散布図を示している中、いくつかの点でatemp=12あたりにtempに拠らず一定値を取る点群がある。
外れ値、というより欠損値が12あたりにアサインされているような雰囲気。
最初の機械学習
最初にシンプルに機械学習モデルLinear Regressionでフィッティングして得られる結果を確認。
結果はKaggleのスコア1.33972。上位92.9%。
工夫無しだとどの程度か？といった確認として。
- 使用モデル：
- scikit-learn提供のLinear Regerssionモデル
- データ分割：
- train用データをKFold交差検証で5分割して「train」と「val」データに分ける
- 結果の評価：
- RMSLE値と、予測値 vs 真値 のプロットを描く
- RMSLE値はscikit-learn提供の「mean_squared_log_error」の平方根より算出する
※scikit-learn ver.1.4以降で「root_mean_squared_log_error」モジュール有り。しかし、2024年5月現在、Kaggle上でのscikit-learnはver.1.2であった為
mean_square_log_errorの日k数オプションsquaredをFalseにすることで値を得る
- 使用データ
- オリジナルのtrain.csvをベースとし、以下変更を加えた変数
- datetimeを年、月、日、時刻、曜日に分けた
- datetime、registered、casual、countを削除した
（途中で作った「Test_Flag」も削除している）
- オリジナルのtrain.csvをベースとし、以下変更を加えた変数
結果
本来直線に乗るべき、予測値 vs 真値 のプロットが直線でなく、$\sqrt{x}$ か log（x）のような分布になっている。
yearの係数が課題。2011年：2012年比は約1.67であったがyearの係数は82。
上手くモデリングできないことが分かる。
ここから改善させていく。
スコア確認1
Kaggleにsubmitしたときのスコア。
- スコア：1.33972
- 順位」3012位 / 3242人
- 上位%：92.9%
ただの線形モデルフィットでも下位200人（下位7%）には勝っているらしい。
下位7%は何かの間違いのような値なのだろうか？
これの逆だとすると、上位7%以上はかなり優秀に思える。
課題
当面思いつく課題。
効率面を考えると遠回りではあるが、まずは要因を書き出して順に潰していくことにする。
- 説明変数の編集
- temp、atempが強い相関を示しており、共線性が問題となる
- atempの分布に歪みがあるので正しいatempを算出する、またはtempのみの使用とする
- 線形モデルなので、seasonや月、時刻等のカテゴリカル変数のエンコード順を昇順または降順にすることで変数の線形近似が良くなると考える
- humidity、windspeedにおいて、他の値から飛んで0となっているものがあり、欠損値を0で埋めたものの可能性がある。0を適正に「欠損値処理」する
- 他、場合により量的データに対して標準化や対数化などの常套手段を検討
- 目的変数の扱い
- casual、registeredで分布の特徴が異なる
- countを直接求めるのではなく、casual、registeredをそれぞれ求めて和を取る
- L2正則化（Ridge回帰）
- 説明変数yearの係数が82と大きい（次に大きい数値は9）。L2正則化を行うRidge回帰が望ましい
説明変数の編集
atempの算出
これは長くなるので別記事にした。
実施内容と得られた結果。
| 実施内容 | 得られた結果 |
|---|---|
| 事実上の欠損値を確認 | humidity=0、atemp=12.12。特定の日付でこの値を示しており、欠損値と判断 |
| ワシントンD.C.の天気予報を取得 | humidity=0を正しく埋める windspeed=0は正しいことを確認 tempのノーマライズ式を仮定、検証 |
| tempのデノーマライズ | 「atemp=12.12」はデノーマライズしたtempの値が0の時のatempの値（切片）であった |
| atemp算出 | デノーマライズしたtemp、windspeed、humidityよりatempを算出 |
| at_norm取得 | atempに施されているノーマライズ式を仮定、検証 算出した正しいatempにもノーマライズを適用 「at_norm」と名付け、欠損無し一部不連続点を排除した正確なatempとした |
カテゴリカル変数のソート
カテゴリカルを線形理論モデルに適用しているのもそもそもどうか、という問題はあるが、
これらをソートすることで線形に準じた形にする。
各カテゴリカル変数を示す。
以下を考慮する。
- 値を3個以上持つカテゴリカル変数に於いてソートが必要
- dayについてはtestデータには値を持たず、予測に使用されない
- weather=4は3つしかエントリーがなく、trainデータに1つ、テストデータに2つのみ
よって、season、month、hour、dayOfWeekについてソートする。
weather=4は誤差程度にしかならない影響度と考え、weather=3に統合した。
hourは先の分析でworkingday=0/1で順序がかなり異なることが分かっているが、workingdayおよびholidayの別の無い状態でソートした。
条件付けすると順序が変わるのは一般的であり他の変数でも同様の為。
df_all_sortというDataFrameを作成し、ソートした（countに対する昇順に振りなおした）変数名を、それぞれの変数名末尾に「_sort」をつけて命名した。
humidity=0, windspeed=0について
これらはatempの算出時に確認できた。
詳細は別記事に記載している。
humidity=0は物理的に有り得ない数字。事実上の欠損値であった。
これは或る1日に於いて（2011年3月10日）のみの欠損であったため、実際の天気予報より取得した。
当該データの実際の値はほぼ100%という値であった。
windspeed=0はほかの値から離れており且つヒストグラムのカウントが大きい為欠損値かと思われたが、正常な値であったと判断する。
実際の天気予報での確認と、一般に風速センサでは最小計測可能速度が定められており、これは風速の分解能に比べて大きな値を示していること、ヒストグラム上でのビン幅を考慮すると不自然ではないと考えられることより結論付けた。
0~6の間が空白になっている。
ほかの変数から簡単に予測値を作ることはできるが（RandomForest等でモデリングすると埋まる）、敢えてここを埋めることは避けた。
微風が目的変数に与える影響はほぼないと考えた。
線形理論でモデリングする際に線形性を多少失うことになるが、予測値で埋めることで本来の値が逆順になることでも同様の影響を与えることになる。
不要に不確定な情報且つそもそも影響度の小さい値を編集することを避けた。
よって、humidity=0は正しく補完を完了、windspeed=0は補完の必要なしとした。
目的変数の扱い
予測対象はcount。
count = registered + casual という関係。
resisteredとcasualとは異なる確率分布を持っているように見られる。
registered、casualをそれぞれ求め、これの和を取ることでcountを算出することとする。
registered、casualの分布は以下。
散布図では、registered vs casualについてworkingdaで色分けするときれいに傾向が分離された。
registeredは勤務日にて高い数値を出し、casualではその逆である。
勤務日か否かが、casualとregisteredの特徴を分かつ項目のようである。
L2正則化（Ridge回帰）
重回帰での係数がyearで82と過大。
L2正則化により係数の肥大化を抑える。
機械学習2
上記説明の「説明変数の編集」、「目的変数の扱い」、「L2正則化」、の3つをそれぞれ試す。
各効果を確認したく、1つずつ適用する。
説明変数の編集
- humidity=0を正しい値に更新
- atemp=12.12(事実上の欠損値)を算出補完
- atempをat_normとして正しい値で算出
- カテゴリカル変数を昇順にソート
以下の不要な情報や、ソート前の変数は削除する
datetime, temp, atemp, t, at, season, month, hour, dayOfWeek, casual, registered
結果の Before - After を記す。
Beforeは最初の機械学習。
散布図は代表例としてKFoldの1番目を表示している。
軸を揃えられていないが容赦願いたい。
| Before | After |
|---|---|
 |
 |
|
 |
|
 |
|
 |
[Good]
スコア$R^2$は約0.39 → 0.65 と改善。
RMSLEも約1.30 → 1.10 と多少改善。
[Bad]
yearの係数が82 → 92とさらに肥大化した。
予測値と真値の散布図は依然$\sqrt{x}$やlog(x)のような形
なお、説明変数の処理、「真の体感温度算出」、「カテゴリカル変数のソート」それぞれ1つずつについて効果を確認した。
| 評価指標 | 元の数値 | 真の体感温度計算 | カテゴリカル変数のソート | 両方実施 |
|---|---|---|---|---|
| RMSLE | 1.302 | 1.301 | 1.097 | 1.096 |
| $R^2$値 | 0.394 | 0.396 | 0.652 | 0.653 |
真の体感温度計算は労力を要したがあまり改善に寄与しない。
カテゴリカル変数のソートは寄与が大きい。
線形理論に於いて順序を揃えると改善するのは当然の結果と理解できる。
カテゴリカル変数を多用することは、線形理論でそもそも適していないのかもしれない。
学習モデルの変更
次に学習モデルの変更を行う。
yearの係数がかなり大きいのでRidge回帰でL2正則化を行う。
ハイパーパラメーターは正則化項にかかる係数、$\alpha$。
0.01～10,000の間で値を試したところ$\alpha = 10000$付近でようやく結果が変わる。
かなり大きい値。
ハイパーパラメーター$\alpha$の最適値はライブラリoptunaを用いてベイズ最適化により取得することにする。
ベイズ最適化の結果
（横軸：$\alpha$の値 縦軸：RMSLE値）
最適値：$\alpha=82081$
結果。Linear Regressionとの比較。
| LinearRegression | Ridge（結果と係数のみ） |
|---|---|
 |
 |
 |
 |
| $R^2:0.651$ | $R^2:0.572$ |
| $RMSLE:1.096$ | $RMSLE:0.971$ |
[Good]
- RMSLE値が1.096→0.971と改善
- yearの係数が90前後から13.6に抑えられた
[Bad]
- $R^2$はLinear Regressionで0.65に対しRidgeで0.572と悪化した
- 結果の散布図を見るとLinear RegressionとRidgeとで特に変化は感じられない
あまり喜べる結果ではなかった。
他、思い立ってカテゴリカル変数のyearを「2011,2012」から「0, 1」にエンコードしたが結果は変わらずであった。
線形モデルでの限界と考える。
そもそも目的変数のヒストグラムを描くと値0が多く右裾長の分布なので最小二乗法で解くモデルは適していないと考えられる。
目的変数の扱い
count = casual + registered であり、casual、registeredはそれぞれ分布が異なることが分かっている。
casual、registeredの値をそれぞれ求めて、最後にこの和からcountを算出することにする。
Ridge回帰でこの効果を確認する。
casual, registered ともにoptunaのベイズ最適化でRidge回帰のハイパーパラメーターを求めたうえで学習を行う。
各分布とスコアを示す。
| casual | registered | count |
|---|---|---|
 |
 |
 |
 |
 |
None |
| $R^2 : 0.497$ | $R^2 : 0.531$ | $R^2 : 0.588$ |
| $RMSLE : 1.058$ | $RMSLE : 0.930$ | $RMSLE : 0.969$ |
[Good]
- RMSLE値が0.971 -> 0.969と改善
- $R^2$が0.572 -> 0.588と改善した
[Bad]
- 結果の散布図を見ると特に変化は感じられない
- casualに於いては$R^2$、RMSLEともにスコアが悪い（0.497、1.058）
casual、registeredに分けることの効果は有る。
やはり線形理論予測困難と思われる。
スコア確認2
Kaggle上でのスコア。
- スコア：1.12338
- 順位」2831位 / 3242人
- 上位%：87.3%
労力をかけた割に全然だめだった。
労力の多くは改善への寄与の低い「atempの正確な値算出」（と、pythonのエラーとの戦い）だったので当然か。
考察と課題
線形モデルでは予測が難しいと考える。
他のモデルを検討。scikit-learnのチートシートを見ると、有力な候補は非線形SVRだが、これはサンプル数（DataFrameの行数）が1万以下で有効。
これ以上の場合計算コストが高くなり（時間がかかり）実用上あまり好まれない。
使用している変数にカテゴリカル変数が多く、求めたいのは期間中の内挿なのでRandom ForestやXGBoost等の決定木での予測が良さそう。
次に、casualで特に予測精度が悪かったことについて。
- ゼロ過剰、右裾長分布：casualで特に顕著。この場合最小二乗法をベースとした予測では誤差を独立で同一の正規分布に従うという仮定の下でフィッティングするためアルゴリズムの前提と合わない。
- RMSLE ： ゼロ付近での多少の誤差がRMSLE値をかなり肥大化させる。casualでは特にゼロが多いためRMSLEでの評価は不利となる。
ゼロを最大とした右裾長の分布では一般化線形モデル、例えばポアソン回帰やゼロ過剰ポアソン回帰、負の二項分布回帰が適用される。
これを使用する場合、アドバンストな内容ではあるが、目的変数の分布に対する検証も必要となる。
線形での予測に限らず、最小二乗法を用いたモデルは厳しいことを念頭に置く必要がある。
LightGBMによる予測
動作の軽い手法として採用。
決定木のブースティングにつき、線形性やカテゴリ変数のソートは気にしなくてよい。
ただし多重共線性の強い説明変数があると説明変数の重要度が不安定になるので、相関の強い説明変数同士はどちらかを削除する。
LightGBMの正しいよう法についてかなり時間をかけて勉強した。
学習内容をベースに考慮した事項を記したいが、多種にわたりここでは書き尽くせない。
LightGBMを理解している者にしか伝わらない内容となるが、肝となる箇所についてのみ記載する。
評価はKFold交差分割検証でデータを5分割、分解したvalidation用データを使用してrmsle評価用の関数を作成してこれで評価した。
同時にKFoldのtrain用データにも同評価を行い、train、validationともに学習曲線を描けるようにした。
KFold交差検証後の本番のモデルでは全データを用い、best_iteration_数だけn_estimatorsを回した。
num_leavesは学習用データの母数に比例する値なので、最適化で得られた値を1.25倍して使用した。
（本番モデル作成ではvalidation用データは確保・使用せず、early_stoppingは使用しない）
LightGBMデフォルト値での学習
最適化前に、一度デフォルト値（パラメータ指定なし）で学習を実施。
パラメータ最適化の結果を見るためのベースとして。
一応、2種類で実施。
- countを直接求めたもの
- casual + registered = count として、casual、registeredを求めたもの
countに対して予測した結果
KFold交差検証
casualに対して予測した結果
KFold交差検証
registeredに対して予測した結果
KFold交差検証
KFold交差検証のRMSLE値平均
| count | casual | registered | |
|---|---|---|---|
| train | 0.372 | 0.494 | 0.365 |
| valid | 0.414 | 0.536 | 0.401 |
実際の予測（trainデータ100%使った予測）
各種RMSLE値はおおよそ0.37～0.5程度であり、LinearRegressionやRidgeの半分以下を示す。
予測値 vs 真値の散布図も直線上に乗ったうえでバラつきを示しているようで、バリアンス大だがバイアスは正しそう。
いい感じ。
スコア確認3
Kaggle上でのスコア。
- スコア：0.51853
- 順位」1878位 / 3242人
- 上位%：57.9%
ようやく半分程度、以前下位。
LightGBM最適化での学習
次にパラメータをチューニングしてLightGBMを使用する。
パラメータはLightGBMの公式ドキュメントと各パラメータの意味より検討した。
なお、LightGBMの実装はsklearn APIをベースにしそのパラメータ名を記載している。
他のAPIに於けるパラメータ名と若干異なるものもあるが、容赦願いたい。
- 最適化対象
- num_leaves
- min_data_in_leaf
- 最適化手法
- optunaによるベイズ最適化
- early_stoppingを使用。n_estimatorsの値はここから得る
- early_stopping
- KFold交差検証でvalidationデータを使用
- validationデータセットに対してrmsleを計算し5ラウンド以上変化しない場合early_stopping
- n_estimatorsの最適値について
- n_estimators はearly_stoppingから得られるbest_iteration_の値を最適値とする
- 評価関数
- モデル学習時のeval_metricには自作したRMSLE値算出関数を設定した
最適化の一例を示す。
optunaの標準関数により、最適値を得るまでの結果をグラフプロットできる。
casualについてのこれを示す。
＜最適化＞
目的変数 = casual、一例として
| 最適化History | 両パラメータ間のコンター図 |
|---|---|
 |
 |
min_data_in_leafが4とかなり小さい。
実際の学習。KFold交差検証。
＜KFold交差検証＞
目的変数 = count
目的変数 = casual
目的変数 = registered
＜最終的な学習結果＞
目的変数 = count
目的変数 = casualとregistered、合計してcountを得る
各RMSLE値
評価のRMSLE値は良くなっている。
- 目的変数=countのとき ： 0.271
- 目的変数=casual、registeredのとき ： 0.246
スコア確認4
Kaggle上でのスコア。
- スコア：0.47609
- 順位」1116位 / 3242人
- 上位%：34.4%
かけた労力はLightGBMの理解のための勉強。
順位は上記3割に迫った。
考察
以下より、過学習と考えられる
- casual、registeredの「min_data_in_leaf」値が4前後であり、かなり小さい
- KFoldのRMLSE値に於いて、trainとvalidの差が大きい
- KFoldの散布図を見ると、trainではデータのバラつきが小さいがvalidでは大きい
- 学習曲線ではvalidが早期に収束しているのに対しtrainでは下降を続けている
最も理解しやすい過学習要因は「min_data_in_leaf=4前後」と極端に小さいこと。
しかし、これはKFoldのvalidationデータを対象に最適化したものなので、これを無暗に大きな値にするのは汎化のみを狙うあまりに視野が狭くなった行為のように思う。
今回のデータを受けて考察され方針出しすべき点は「汎化に向けて」という結論になると思われるが、一度これは保留とする。
実はmin_data_in_leafが小さくなるのはモデルの設定ミスや説明変数の問題等ではなく、
目的変数の分布がゼロで最大を示し右裾長になっているので仕方がないと考えることもできる。
数学的に説明するほど理解できていないのだが、滅多に起こらない事象が多くを占める分布に対して決定木を作ろうとすると、かなり分岐しないとレアな事象を特定できないと思われる。
他に考えられる施策を行う。
スコアを上げるための施策
スコアを上げるための方針を2点実施する。
- 説明変数の作成
- 目的変数の分布考察
説明変数の作成
予測するための情報ついて考える。
予測するうえで、通常は「ユーザ情報」が重要になる。
よく見るユーザアンケートでは年齢、性別、利用目的等のユーザの属性を基に傾向を予測することが多い。
今回の説明変数では気象条件と日時のみの用意となる。
また、Kaggleコンペ色が強い考え方だが、予測結果はRMLSEで評価されるので小さな値を示す箇所で正確に予想するとスコアが伸びやすい。
よって、
- ユーザ属性を反映した説明変数を作る
- 利用者数の低い個所の特徴を見る
ということを行いたい。
通勤利用者、レジャー利用者
workingday=1とworkingday=1とで分布が異なる。
registeredの分布をみると通勤利用者の像がイメージできる。
これを基に特徴量を作成する。
| 特徴名 | 値 | 条件 |
|---|---|---|
| commute | 0,1 | workingday = 1 , hour = 7 or 8 or 17 or 18 |
| leisure | 0,1 | workingday = 0 , hour = 11 to 17 |
付随して、レジャー利用者というものも付け足した。
本当は「祝日でも働いており、通勤利用している人」がいるように見える。
ユーザ属性で、「祝日関係ない業務に従事している」というユーザ層が分かれば精度が上がりそうなのだが、残念。
1日雨が続く
本コンペのデータは、時系列データとしての性格を強く持っている。
一度時系列でデータを見てみたい。
データ量が多くかつ時刻に依って乱高下するので、日単位でグループ平均してプロットした。
workingday、holiday、weatherも併せてプロットする。
これらは上下動が視認できるように、適当にスケーリングして表示している。
ところどころ急峻な落ち込みを示す箇所がある。
日でグループ平均したweatherと関係が強い。
weatherは瞬時値を示しておりこれとは正確に対応していないようなので、
長期間にわたる悪天候を特徴化したい。
「今日は一日中雨」という長期の系と、「昼だけ雨」という瞬時値の系とではレンタルバイク利用に与える影響度は異なりそう。
グラフ作成に倣い、precipitation（降雨量）としてweather値の24時間移動平均を作成する（重み無し平均）。
現在時刻に対する荷重が高いかと思い重み付けを試みたが、結果的にローパスフィルタの効果が急峻になるだけなので、重み無しの移動平均とした。
precipitation > 2 となるデータにbad_weather という0/1の値を付与する。
日単位グループ平均でbad_weather=1に該当する個所を赤で塗りつぶして表示する。
（無論、日（24時間）単位で平均された表示となるのでweatherの値とbad_weatherとの値はほぼ変わらないが、効果を俯瞰から確認するという意味で理解してほしい）
急峻な落ち込みの多くをカバーしているように思う。
時刻単位の例は、代表2日を表示。
bad_weather=1に該当するエリアを赤色で塗りつぶして表示した。
休祝日ごとに落ち込みがあるので分かりにくいが、うまく特徴をとらえていると思う。
また、bad_weather（およびprecipitation）にかからない落ち込みはweather=3の瞬時値が捉えている。
作成した特徴量まとめ。
| 特徴名 | 値 | 条件 |
|---|---|---|
| bad_weather | 0,1 | precipitation>2 |
*precipitation : weatherの24時間移動平均
利用者数が多い条件
precipitatoin、bad_weatherで利用者数が極小になる条件を特徴化した。
折角なので、これの逆もやってみたい。
precipitation作成時に作図した日単位グループの時系列データにて、
極大ピークを示す点（年月日）を15個ピックアップ。
その条件を見て共通条件（条件範囲）を調べた。
特筆すべきはworkingday=0が多く特にcasualの値が大きいという条件が目立った。
これをcomfortと命名し、以下に特徴化する。
| 特徴名 | 値 | 条件 |
|---|---|---|
| comfort | 0,1 | 18 < at_norm < 35 & 40 < humidity < 80 & windspeed < 18 & precipitation < 1.7 & workingday = 0 |
クリスマス
よく見ると、クリスマスがholidayになっていなかった。
（厳密にいうとクリスマス前後で休暇になっていたり、なっていなかったり。
2011年と2012年とで扱いが異なるようであった。）
欧米では年始の正月休みがない代わりにクリスマス以降休暇に入る文化があると経験している。
ドイツ系企業で勤務していたが、12/23を最終日として、以降連絡がつかなくなる。
アメリカでも同様であろう。
思い切って12月24日以降はholiday=1にしてしまう。
workingdayも併せて0にする。
monthを削除
日単位グループ平均のデータを見ると、特にbad_weatherの関係で落ち込む月がちらほらある。
本記事の序盤で、「monthは2011年、2012年とで傾向は同じ」と述べたが、よく見ると多少異なる。
月情報は19日間分の訓練用データ情報で、うち5日間程度悪天候が続く月などがある。
2011年と2012年とでこの様子は同様ではないことに気が付く。
以下、monthの分布。各年のcount量に対する各月の占める割合を示す。
seasonでは以下の通り。
monthとseasonとは相関の強い2変数だが、monthを削除してseasonを採用し、情報の平準化をしたほうが良いと判断した。
作成した説明変数まとめ
| 変数名 | 表現する特徴 | 条件設定 |
|---|---|---|
| commute | 通勤利用者 | workingday=1の通勤ピーク時刻 |
| leisure | レジャー利用者 | workingday=0の日中 |
| precipitation | 利用者が極小となる条件 | 長時間悪天候を示す |
| comfort | 極大となる箇所、casualをベースに作成 | ピーク個所の気象条件を並べた |
| その他 | クリスマスを休暇に | 12月24日以降をholidayに設定 |
＝＝＝＝＝以下、その具体的設定値を再度記載＝＝＝＝＝
| 特徴名 | 値 | 条件 |
|---|---|---|
| commute | 0,1 | workingday = 1 , hour = 7 or 8 or 17 or 18 |
| leisure | 0,1 | workingday = 0 , hour = 11 to 17 |
| 特徴名 | 値 | 条件 |
|---|---|---|
| bad_weather | 0,1 | precipitation>2 |
*precipitation : weatherの24時間移動平均
| 特徴名 | 値 | 条件 |
|---|---|---|
| comfort | 0,1 | 18 < at_norm < 35 & 40 < humidity < 80 & windspeed < 18 & precipitation < 1.7 & workingday = 0 |
その他、クリスマスを休暇に設定。workingdayも併せて0に。
casual、registered での使用
| casualで使用 | registeredで使用 |
|---|---|
| precipitation | precipitation |
| --- | commute |
| leisure | leisure |
| bad_weather | --- |
| comfort | --- |
それぞれ、casual、registeredの特徴を記したものはもう一方の特徴量にあてはめが効かないようで、
実験的に交差検証したところ上記組み合わせで使用するとスコアが伸びた。
もともとの変数作成モチベーションと照らし合わせると合点がいく結果である。
スコア確認5
Kaggle上でのスコア。
スコア：0.42761
順位」536位 / 3242人
上位%：16.5%
順位は前回の半分になった。
特徴量の追加は非常によく働いた。
目的変数の分布について
目的変数の分布について。何度か触れているが再度。
特にcasualのヒストグラムで顕著だが、0が最大となり右に裾が長い分布となる。
この場合、正規分布を仮定する重回帰分析の前提条件から外れる。
（最小二乗法の観点から、正確には撹乱項(残差の確率分布項)が独立に同一の正規分布に従うことを仮定）
casualのヒストグラム
registeredのヒストグラム
registeredは、workingday=0のときcasualに類似
一方で、ポアソン分布に類似した形を示しているが、分布形状を比べると異なる。
ポアソン分布の特徴は**平均値と分散とが同じ値（$\lambda$）**となる。
casualについてはdescribe()関数で手軽に平均meanと標準偏差stdが得られる。
casualでは、平均（mean）と分散（std値を2乗した値）は大きく異なっており、ポアソン分布には従わないことが分かる。
（※メモ：DataFrameのdescribe()関数は標本標準偏差（n-1で割る）なので、ここで議論したい値と比べると若干値は大きめに記される。
一様分布であるdayOfWeek等に於いてmeanとstdとの関係が合わないのはこのためと思われる）
最尤度推定により特別な確率分布を求める一般化線形モデルという手法もある。
ポアソン回帰やゼロ過剰XXX、負の二項分布がその代表例。
これは推定に使う確率分布の妥当性を検証することが肝となるが、今回は分布が難しい。
先に作成した「commute」の条件でヒストグラムを作成する。
通勤利用者に絞ってフィルターしたものである（平日の7時、8時、17時、18時に於けるregisteredのユーザ）。
registeredのhourに於けるbarplotの最大値を示すエリアであったので、registered全体のヒストグラムの0を除いたピークである150前後の分布を形作るのかと思ったが、450あたりに最大を持つ分布であった。
また、ピークらしき点が3か所見られる。
気象条件や季節によりピークが前後する、または異なるユーザ区分が混ざっておりそれぞれの分布が合成される、など考えられる。
以下はcommute=1（registeredのピーク条件を抽出）に於ける時刻歴波形。
経時的なゆるやかな変化と、急峻な落ち込みがみられる。
天候が悪い条件の日時に赤塗、緑塗りを施すと多くはこれに該当する模様。
一般に、ポアソン回帰等行う場合には急峻な落ち込みは外れ値処理を行い分布を整えて実行する。
今回は外れ値というには条件が定まっており且つ頻度も高いため予測に使用したく、外れ値処理にかけたくない。
ポアソン回帰にゼロ過剰を考慮したものや負の二項分布での回帰等、確率分布を利用した回帰手法の検討をしたが、上記の通り確率分布の調査の段階でこれを棄却した。
副次的に、一般化回帰に使用するリンク関数として、目的変数にlogを取る手法がある。
今回はこれを採用することにする。
ただし、目的変数に0が多数あるため、真数条件を考慮しlog(x+1)で変換する。
これのヒストグラムを描くと分布は極端な右裾長ではなくなる。
この分布を以て機械学習を行うと比較的多くのエリアで自然な分布（正規分布）に近い状況に対する予測に落とし込むことができる。
予測能力向上を大きく見込めそうである。
スコア確認6
Kaggle上でのスコア。
スコア：0.37525
順位」94位 / 3242人
上位%：2.9%
飛躍的に改善。上位3%以上になった。
一旦これで目標は達成したことになる。
この時のLightGBMのパラメータを確認すると、以下の通り。
min_data_in_leafは最初の最適化時には2や4という極小さい値であり過学習が懸念された。
最終的には微量ながらこれも改善されており、casualで6、registeredで23となった。
LightGBM パラメータ
＜casual＞
| ハイパーパラメータ | 値 |
|---|---|
| n_estimators | 74 |
| num_leaves | 130 |
| min_data_in_leaf | 6 |
＜registered＞
| ハイパーパラメータ | 値 |
|---|---|
| n_estimators | 107 |
| num_leaves | 128 |
| min_data_in_leaf | 23 |
まとめ
KaggleのBike Sharing Demand に予測を行った。
本コンペの内容の特徴は以下の通り。
- 予測する期間は内挿：毎月20日以降の予測
- 評価値はRMSLE値
- 目的変数の分布は0が多い右裾長の分布
- 目的変数はcasual、registeredという2種類の属性の和となる
- 説明変数は気象条件と年月日、平日、休日情報のみで利用者に関する情報はない
- 見かけ上欠損値は無いが、事実上の欠損値がある
分析時の特徴は以下の通り。
- casual、registeredで特色の異なる分布を示す
- 目的変数をlog(x+1)とすると機械学習の精度を上げやすい
- 決定木をベースとしたLightGBMで比較的良い予測、線形モデルでは分布の形状に追従できない（バリアンスがかなり高くなる）
- 利用者に関する情報がない為、分布からそれに相当するものを作ると良い
考察
分布がゼロ過剰なうえ、評価指標がRMSLEなのでスコアを上げることが難しい問題であった。
Kaggleリーダーボード上の1位のスコアも0.33756であり、凡その目安でプラス側40%、マイナス側30%の誤差量に相当。
一見すると悪いスコアだが、RMSLEではゼロ付近の誤差が大きく響く特徴を持ち、今回の分布はゼロが多いことに依ると考えられる。
今回の自己のスコアは0.37525、上位約2.9%。
「機械学習」を学んで3か月程度で取得した順位としては悪くないと思う。
一方、モデルは簡易にLightGBMを使用し特徴量と目的変数について考慮したのみ。
データに関する解釈性については現状で良いが、さらに一歩踏み込んでスコアアップをすることもできそう。
アンサンブルにより弱学習機の結果を合わせることでバリアンスを抑えることを試してみたい。
中心極限定理（平均値の寄り集め集団と母平均との誤差は正規分布に収束する）のように良い値に収束していくのかと思われる。
今回のモデルでは予測値と真値の散布図をプロットするとバイアスはあまり感じられないようなので、そこそこ良い結果を得られそう。
追記
→後日、アンサンブルを行った。
これまでのEDAより線形モデルでもR2スコア0.8程度の量スコアを出すようになっており係数も$\pm2$以内に収まっていた。
Linear Regression、Linear SVR、SVR、Random Forest、LightGBMの5個でアンサンブルを行ったがスコアは悪くなった。
唯一、Random ForestとLightGBMの2つのみを2:8で平均するとRMSLE=3.7422、上位2.84%になった。
感触的に、恐らくモデルを20個ほど用意して初めて効果が出てきそう。
そこまでの多様なモデルは容易に作れなかった。
最後に、Kaggleコンペと実務について。
コンペではスコアを競うという特性であり、今回の問題設定は現実のタスクとは離れた問題設定となっている。
各月の20日以降の予測という問題設定だが、実際のタスクでは未来の予測をしたいと思われる。
この場合、2011年→2012年で利用者数が約1.67倍増加していることより、2013年の予測もかなり利用者数は増加しつつ気象や休祝日の影響を受けて変動するものと考えられる。
LightGBMによる決定木ベースのモデルでは外挿には弱いと考えており、また、利用者数増加量については今回のモデルでは「year」特徴量がその解釈を担っていると思われる。
未来を予測したい場合、year相当の特徴量は年間の統計がすべて出た後に定義可能なので使用できない。
時系列データでの分析から増加量の傾向を特徴化して、これに気象や休祝日による変動の効果を考慮して、といったモデルを作成できると良いように思われる。
こうした分析も検討したい。
今回は記事が冗長になること、Aidemy様講座の最終課題として目標とした結果を提出することよりここで終了とする。
後日上記の追加モデルを検討したい。
その前に、離職中なので就職活動に勤しまなければいけない。
雇ってくれる方募集。