目次

自己紹介
実行環境
成果物の概要
ライブラリのインポート
データの読み込み
データの確認
前処理と特徴量エンジニアリング
反省と今後

自己紹介

現在、主婦ときどきフリーランスWEBデザイナーとして細々と仕事をしています。

元々は文系で営業職だった私が、結婚を機に離職し、主人の勧めもあってIllustratorが使えるようになりたくてWEBデザインの職業訓練に通い始めたのがキャリアチェンジのきっかけでした。
ありがたいことにパートタイマーとしてEC関連企業で2年ほどお仕事させてもらいましたが、成長に限界を感じ離職、スキルアップ期間をとることにしました。

デザインというよりもHTML等のコーディングのほうがおもしろかったこともあり、プログラミングに興味を持ちました。
そこから、2023年9月からAidemyのデータ分析講座(6か月コース)を受講するに至ります。

おかげさまで受講3か月ほどで成果物作成の段階に辿り着き、ここに公開します。
なお、このブログはAidemy Premiumのカリキュラムの一環で、受講修了条件を満たすために公開しています。

実行環境

Google Colaboratory

成果物の概要

Store Sales – Time Series Forecasting
https://www.kaggle.com/competitions/store-sales-time-series-forecasting/overview

コンペの目的

時系列予測を使用して、エクアドルを拠点とする大手食料品小売業者 Corporación Favorita のデータに基づいて店舗の売上を予測します。

評価基準

コンペの評価指標は、RMSLE（二乗平均平方根対数誤差）です。

RMSLE=\sqrt{ \frac{1}{n} \sum_{i=1}^n \left(\log (1 + \hat{y}_i) - \log (1 + y_i) \right)^2}

$n$ は、インスタンスの合計数です。
$\hat{y}_i$は、インスタンス (i) のターゲットの予測値です。
$y_i$は、インスタンス (i) のターゲットの実際の値、および、
$\log$は自然対数です。

提出ファイルフォーマット

テストセット内のIDごとに、 sales変数の値を予測する必要があります。ファイルは次の形式である必要があります。

id,sales
3000888,0.0
3000889,0.0
3000890,0.0
etc.

データについて

train.csv

トレーニングデータ。時系列の特徴store_nbr、 family、およびonpromotionと目標salesで構成されています。

store_nbr は、製品が販売されている店舗を識別します。
familyは販売される製品の種類を識別します。
sales は、特定の日付における特定の店舗の製品ファミリーの総売上高を示します。
onpromotion は、特定の日付に店舗でプロモーションされていた製品ファミリー内のアイテムの総数を示します。

test.csv

トレーニングデータと同じ機能を持つテストデータ。このファイルで、日付の目標売上を予測します。テストデータの日付は、トレーニングデータの最後の日付から15日間です。

stores.csv

city、state、type、clusterなどのメタデータを保存します。clusterとは、類似したストアをグループ化したものです。

【事前案内】
上記以外に、以下のようなデータも用意されています。
本来であればこのようなデータも活用し予測を行うべきですが、今回はtrainデータから一部の期間を抽出したデータとtestデータ、またstoresのデータを用いてシンプルな予測をします。

transactions.csv

日付、store_nbr、取引量で構成されています。

oil.csv

毎日の石油価格。トレーニングデータとテストデータの両方の時間枠での値が含まれます。

holidays_events.csv

休日とイベントのメタデータ。

ライブラリのインポート

import numpy as np
import pandas as pd

from sklearn.linear_model import LinearRegression

import warnings

データの読み込み

train = pd.read_csv("/content/train.csv")
test = pd.read_csv("/content/test.csv")
stores = pd.read_csv("/content/stores.csv")

データの確認

一つずつ、ざっくりどのようなデータなのか、確認していきます。

まず、確認しやすい表示のするため、下記のような関数を準備します。

def basic_eda(df):
    print("\n---Head 5---\n")
    print(df.head(5))
    print("\n---Columns---\n")
    print(df.columns)
    print("\n---NULL---\n")
    print(df.isna().sum())
    print("\n---Shape---\n")
    print(df.shape)
    numeric_variables = df.select_dtypes(include=[np.number])
    categorical_variables = df.select_dtypes(exclude=[np.number])

    print('Numeric variables in the given the dataframe : ',numeric_variables.shape[1])
    print('Categorical variables in the givne the dataframe:',categorical_variables.shape[1])

作成した関数を使って、それぞれのデータを見ていきます。

print("\n==== train ====\n")
basic_eda(train)
print("\n==== test ====\n")
basic_eda(test)
print("\n==== Shops ====\n")
basic_eda(stores)

trainデータは2013年からのデータがあり、testデータは2017年8月16日からで、上に説明のあったとおりこの日から15日間の売上(sales)を予測することとなります。

train.date = pd.to_datetime(train.date)
test.date = pd.to_datetime(test.date)

日付のデータ型がobjectになっていたので、datetime型に変換します。

前処理と特徴量エンジニアリング

まずは特徴量を選び、featuresに格納します。目的変数はsalesです。

ここに、日付も特徴量として追加したいので、曜日・月・年を抽出してfeaturesに追加していきます。

features = ['store_nbr', 'family', 'onpromotion']
target = 'sales'

def extract_weekday(df):
    return df.dayofweek

def extract_month(df):
    return df.month

def extract_year(df):
    return df.year

train['weekday'] = train['date'].apply(extract_weekday)
train['year'] = train['date'].apply(extract_year)
train['month'] = train['date'].apply(extract_month)

test['weekday'] = test['date'].apply(extract_weekday)
test['year'] = test['date'].apply(extract_year)
test['month'] = test['date'].apply(extract_month)

features.append('weekday')
features.append('year')
features.append('month')

訓練データとテストデータを、それぞれ'store_nbr', 'family', 'date', 'sales'の列を使用し、データ型の指定を行い、データを整えます。整えた後のデータは、それぞれstore_salesとtest_dfにしました。

store_sales = pd.read_csv(
    '/content/train.csv',
    usecols=['store_nbr', 'family', 'date', 'sales'],
    dtype={
        'store_nbr': 'category',
        'family': 'category',
        'sales': 'float32',
    },
    parse_dates=['date'],
    infer_datetime_format=True,
)
store_sales['date'] = store_sales.date.dt.to_period('D')
store_sales = store_sales.set_index(['store_nbr', 'family', 'date']).sort_index()


test_df = pd.read_csv(
    '/content/test.csv',
    dtype={
        'store_nbr': 'category',
        'family': 'category',
        'onpromotion': 'uint32',
    },
    parse_dates=['date'],
    infer_datetime_format=True,
)
test_df['date'] = test_df.date.dt.to_period('D')
test_df = test_df.set_index(['store_nbr', 'family', 'date']).sort_index()

今回はシンプルな予測になるので、
テストデータは2017年8月16日からなので、訓練データの目的変数は2017年のデータに絞ります。
なお、先ほど作成したstore_salesは下記のようなデータになっています。

display(store_sales)

ここから、日別に売上が見れるようunstackし、2017年分を抽出します。

# Target
y = store_sales.unstack(['store_nbr', 'family']).loc["2017"]

訓練データを作成します。

モデリング（線形回帰）

今回は、線形回帰でどこまでの制度が出るか見てみます。

model = LinearRegression(fit_intercept=False)
model.fit(X, y)
y_pred = pd.DataFrame(model.predict(X), index=X.index, columns=y.columns)

X_test = dp.out_of_sample(steps=16)
X_test.index.name = 'date'
X_test['NewYear'] = (X_test.index.dayofyear == 1)

RMSEによるモデルの精度評価を行います。

冒頭でも述べたように、RMSLEでの制度評価を行うためには、実際の予測値（y_pred）と真の値（y）の対数を取って、その差を評価する必要があります。

y_pred_clipped = np.clip(y_pred, 1e-10, None)
y_clipped = np.clip(y, 1e-10, None)

y_log = np.log1p(y_clipped)
y_pred_log = np.log1p(y_pred_clipped)

log_diff = y_log - y_pred_log

log_diff = np.where(log_diff < 0, 0, log_diff)

rmsle = np.sqrt(np.mean(log_diff**2))

print(f'RMSLE: {rmsle}')

反省と今後

今回の作業で、一通りの流れを体験することができました。
この作業を通して初めて、受講当時にいろんな人から聞いていた、「データサイエンスはほとんどが前処理で地道な作業」という話に納得することができました。
実際、このブログに着手するとなって最初にしていたことはデータとひたすらにらめっこする時間でした。

結果、今回はほとんどのデータを活用することができていません。
今後は受講修了後もこの内容に引き続きチャレンジして、記事をブラッシュアップしていきたいなと思います。