{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Favorita - 重回帰分析\n\n**モデル**: 最小二乗法（OLS）重回帰分析  \n$$\\hat{y} = \\beta_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\cdots + \\beta_n x_n$$","metadata":{}},{"cell_type":"code","source":"import os, subprocess, gc, warnings, time\nimport pandas as pd\nimport numpy as np\nfrom sklearn.linear_model import LinearRegression\nwarnings.filterwarnings('ignore')\n\nINPUT_DIR = '/kaggle/input/competitions/favorita-grocery-sales-forecasting'\nWORK_DIR  = '/kaggle/working'\n\nt0 = time.time()\ndef elapsed(): return f\"{time.time()-t0:.0f}s\"\n\n# 7z 展開\ndef find_7z(name):\n    \"\"\"INPUT_DIR以下から {name}.7z を探す\"\"\"\n    for root, _, files in os.walk(INPUT_DIR):\n        for f in files:\n            if f == f'{name}.7z':\n                return os.path.join(root, f)\n    raise FileNotFoundError(f'{name}.7z が見つかりません: {INPUT_DIR}')\n\ndef extract_7z(src, dst_dir=WORK_DIR):\n    \"\"\"7z ファイルを dst_dir に展開 (7z CLI → py7zr の順で試行)\"\"\"\n    os.makedirs(dst_dir, exist_ok=True)\n    for cmd in [['7z','e'], ['7za','e']]:\n        try:\n            r = subprocess.run(cmd + [src, f'-o{dst_dir}', '-y'],\n                               capture_output=True, timeout=600)\n            if r.returncode == 0:\n                return\n        except (FileNotFoundError, subprocess.TimeoutExpired):\n            pass\n    import py7zr\n    with py7zr.SevenZipFile(src, 'r') as z:\n        z.extractall(path=dst_dir)\n\n# 小さファイルを一括展開\nSMALL_FILES = ['items.csv', 'stores.csv', 'oil.csv',\n               'holidays_events.csv', 'test.csv']\n\nprint('ファイルを展開中...')\nfor name in SMALL_FILES:\n    dst = os.path.join(WORK_DIR, name)\n    if not os.path.exists(dst):\n        src = find_7z(name)\n        extract_7z(src)\n        print(f'  {name} 展開完了 ({os.path.getsize(dst)/1e6:.1f} MB)')\n    else:\n        print(f'  {name} 既存')\n\nprint(f'完了 [{elapsed()}]')\n\n# 補助テーブル読み込み\nW = WORK_DIR\ntest     = pd.read_csv(f'{W}/test.csv',              parse_dates=['date'])\nitems    = pd.read_csv(f'{W}/items.csv')\nstores   = pd.read_csv(f'{W}/stores.csv')\noil      = pd.read_csv(f'{W}/oil.csv',               parse_dates=['date'])\nholidays = pd.read_csv(f'{W}/holidays_events.csv',   parse_dates=['date'])\n\nprint(f'test: {test.shape}  items: {items.shape}  stores: {stores.shape}  [{elapsed()}]')\n\n# train.csv 展開 → チャンク読み込み → 展開ファイル削除\nCUTOFF = '2017-05-16'  # 一旦　lag_16 のみ学習 (学習開始: 2017-06-01)\ntrain_csv = os.path.join(WORK_DIR, 'train.csv')\n\nif not os.path.exists(train_csv):\n    print('train.csv.7z を展開中 (数分かかります)...')\n    extract_7z(find_7z('train.csv'))\n    print(f'展開完了: {os.path.getsize(train_csv)/1e9:.2f} GB  [{elapsed()}]')\n\nprint('train.csv をチャンク読み込み中...')\nchunks = []\nfor chunk in pd.read_csv(train_csv, parse_dates=['date'],\n                         dtype={'unit_sales':'float32','onpromotion':'str'},\n                         chunksize=10_000_000, low_memory=False):\n    chunk = chunk[chunk['date'] >= CUTOFF]\n    if not len(chunk): continue\n    chunk['unit_sales'] = chunk['unit_sales'].clip(0)\n    chunk['onpromotion'] = chunk['onpromotion'].map({'True':1,'False':0}).fillna(0).astype('int8')\n    chunks.append(chunk)\ntrain = pd.concat(chunks, ignore_index=True)\ndel chunks; gc.collect()\n\n# ディスク節約のため展開済みの train.csv を削除\nos.remove(train_csv)\nprint(f'train.csv 削除 (ディスク節約)')\nprint(f'train: {train.shape}  [{elapsed()}]')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T07:41:39.012702Z","iopub.execute_input":"2026-06-01T07:41:39.013038Z","iopub.status.idle":"2026-06-01T07:44:39.144762Z","shell.execute_reply.started":"2026-06-01T07:41:39.013002Z","shell.execute_reply":"2026-06-01T07:44:39.143814Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 補助テーブル前処理\n\n# 原油価格: 線形補間\noil_r = pd.DataFrame({'date': pd.date_range(oil['date'].min(), '2017-08-31')})\noil_f = oil_r.merge(oil, on='date', how='left')\noil_f['dcoilwtico'] = oil_f['dcoilwtico'].interpolate('linear').astype('float32')\ndel oil, oil_r\n\n# 国祝日フラグ\nhol_dates = set(\n    holidays[(holidays['locale']=='National') & (~holidays['transferred'])]['date']\n)\ndel holidays\n\n# 店舗タイプ\nstores['type_enc'] = stores['type'].map({'A':0,'B':1,'C':2,'D':3,'E':4}).astype('int8')\n\n# 商品カテゴリ\nfam_map = {f:i for i,f in enumerate(sorted(items['family'].unique()))}\nitems['family_enc'] = items['family'].map(fam_map).astype('int8')\n\n# 日付特徴量\ndef add_date_features(df):\n    df = df.copy()\n    df['dow']        = df['date'].dt.dayofweek.astype('int8')   # 0=月 〜 6=日\n    df['month']      = df['date'].dt.month.astype('int8')\n    df['is_weekend'] = (df['dow'] >= 5).astype('int8')\n    df['is_holiday'] = df['date'].isin(hol_dates).astype('int8')\n    df = df.merge(oil_f[['date','dcoilwtico']], on='date', how='left')\n    df = df.merge(stores[['store_nbr','type_enc','cluster']], on='store_nbr', how='left')\n    df = df.merge(items[['item_nbr','family_enc','perishable']], on='item_nbr', how='left')\n    return df\n\ntrain = add_date_features(train)\ntest  = add_date_features(test)\n\n# ラグ特徴量: lag_16 のみ (merge-based)\n# 「16日前の売上」を特徴量として追加\nsales_hist = train[['date','store_nbr','item_nbr','unit_sales']].copy()\nsales_hist['date'] = sales_hist['date'] + pd.Timedelta(16, 'D')  # 16日ずらす\nsales_hist.rename(columns={'unit_sales':'lag_16'}, inplace=True)\n\ntrain = train.merge(sales_hist, on=['date','store_nbr','item_nbr'], how='left')\ntest  = test.merge(sales_hist,  on=['date','store_nbr','item_nbr'], how='left')\ndel sales_hist; gc.collect()\n\nprint(f'特徴量作成完了  train: {train.shape}  [{elapsed()}]')\nprint(train[['date','dow','month','is_holiday','dcoilwtico','type_enc',\n             'family_enc','perishable','onpromotion','lag_16','unit_sales']].head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T07:44:39.147010Z","iopub.execute_input":"2026-06-01T07:44:39.147307Z","iopub.status.idle":"2026-06-01T07:44:51.413210Z","shell.execute_reply.started":"2026-06-01T07:44:39.147281Z","shell.execute_reply":"2026-06-01T07:44:51.412310Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"FEATURES = [\n    'lag_16',       # 16日前売上\n    'dow',          # 曜日\n    'month',        # 月\n    'is_weekend',   # 週末フラグ\n    'is_holiday',   # 祝日フラグ\n    'onpromotion',  # 販促フラグ\n    'dcoilwtico',   # 原油価格\n    'type_enc',     # 店舗タイプ\n    'cluster',      # 店舗クラスタ\n    'family_enc',   # 商品カテゴリ\n    'perishable',   # 生鮮品フラグ\n]\n\n# 学習期間: 2017-06-01 以降 (直近2ヶ月)\nTRAIN_START = '2017-06-01'\nVAL_START   = '2017-07-31'\n\nmask_tr  = (train['date'] >= TRAIN_START) & (train['date'] < VAL_START)\nmask_val = (train['date'] >= VAL_START)\n\n# NaN を 0 で補完 (lag が存在しない行)\nX_tr  = train.loc[mask_tr,  FEATURES].fillna(0).values.astype('float32')\ny_tr  = np.log1p(train.loc[mask_tr,  'unit_sales'].values)   # log変換\nX_val = train.loc[mask_val, FEATURES].fillna(0).values.astype('float32')\ny_val = np.log1p(train.loc[mask_val, 'unit_sales'].values)\nw_val = np.where(train.loc[mask_val, 'perishable'].values == 1, 1.25, 1.0)\n\nprint(f'学習: {X_tr.shape}  バリデ: {X_val.shape}  [{elapsed()}]')\n\n# OLS 重回帰分析\nlr = LinearRegression(n_jobs=-1)\nlr.fit(X_tr, y_tr)\nprint(f'学習完了  [{elapsed()}]')\n\n# 回帰係数の表示\ncoef_df = pd.DataFrame({\n    '特徴量': FEATURES,\n    '係数':   lr.coef_\n}).sort_values('係数', ascending=False)\n\nprint(f'\\n切片 β₀ = {lr.intercept_:.4f}')\nprint('\\n回帰係数:')\nprint(coef_df.to_string(index=False))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T07:44:51.414260Z","iopub.execute_input":"2026-06-01T07:44:51.414544Z","iopub.status.idle":"2026-06-01T07:44:53.683039Z","shell.execute_reply.started":"2026-06-01T07:44:51.414519Z","shell.execute_reply":"2026-06-01T07:44:53.681705Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# バリデーションスコア\n\ndef nwrmsle(y_true_log, y_pred_log, weights):\n    diff = y_pred_log - y_true_log\n    return float(np.sqrt(np.sum(weights * diff**2) / np.sum(weights)))\n\nval_pred = lr.predict(X_val).clip(0)\nscore = nwrmsle(y_val, val_pred, w_val)\nprint(f'Local NWRMSLE (OLS重回帰): {score:.5f}')\n\n# R² (決定係数)\nss_res = np.sum((y_val - val_pred) ** 2)\nss_tot = np.sum((y_val - y_val.mean()) ** 2)\nr2 = 1 - ss_res / ss_tot\nprint(f'R² (決定係数):              {r2:.4f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T07:50:02.432186Z","iopub.execute_input":"2026-06-01T07:50:02.432589Z","iopub.status.idle":"2026-06-01T07:50:02.492773Z","shell.execute_reply.started":"2026-06-01T07:50:02.432562Z","shell.execute_reply":"2026-06-01T07:50:02.491856Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_te = test[FEATURES].fillna(0).values.astype('float32')\n\ntest_pred = np.expm1(lr.predict(X_te)).clip(0, 1000) # log変換を戻す & 負値除去\n\nsub = test[['id']].copy()\nsub['unit_sales'] = test_pred\n\nout_path = f'{WORK_DIR}/submission_ols.csv'\nsub.to_csv(out_path, index=False)\n\nprint(f'submission_ols.csv 保存完了: {out_path}')\nprint(f'行数: {len(sub)}')\nprint(sub['unit_sales'].describe())\nprint(f'\\n総経過時間: {elapsed()}')\n\nprint(f'''\n── 結果まとめ ────────────────────────────────────\nモデル        : OLS 重回帰分析\nLocal NWRMSLE : {score:.5f}\nR²            : {r2:.4f}\n学習期間      : {TRAIN_START} 〜 {VAL_START} 直前\n特徴量数      : {len(FEATURES)}\n\n───────────────────────────────────────────────────\n''')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-01T07:44:54.715556Z","iopub.execute_input":"2026-06-01T07:44:54.715834Z","iopub.status.idle":"2026-06-01T07:45:03.263543Z","shell.execute_reply.started":"2026-06-01T07:44:54.715808Z","shell.execute_reply":"2026-06-01T07:45:03.262689Z"}},"outputs":[],"execution_count":null}]}