{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":9994,"databundleVersionId":752467,"sourceType":"competition"}],"dockerImageVersionId":31259,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:31:51.448171Z","iopub.execute_input":"2026-01-20T08:31:51.448507Z","iopub.status.idle":"2026-01-20T08:31:51.470616Z","shell.execute_reply.started":"2026-01-20T08:31:51.448480Z","shell.execute_reply":"2026-01-20T08:31:51.469780Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Monte Carlo Uncertainty Propagation with LMS for Building Energy Modeling\n\n**From deterministic predictions to risk-aware energy forecasting**\n\nReal-world energy forecasting systems operate on noisy, imperfect measurements.\nTreating these inputs as exact values leads to overconfident models that fail in production.\n\nThis notebook demonstrates how **Monte Carlo simulation can be layered on top of an online LMS (SGD) regression model**\nto quantify **prediction uncertainty** for building-level energy consumption.\n\n> **Key idea:** Two predictions with the same value can have very different levels of confidence.\n\nGitHub repository:\nhttps://github.com/HarshDubey-12/monte-carlo-lms-energy-modeling\n","metadata":{}},{"cell_type":"markdown","source":"## Dataset Overview\n\nWe use building-level energy data from the ASHRAE dataset:\n- **Energy meter readings** (target variable)\n- **Building metadata** (static context)\n- **Weather data** (uncertain exogenous inputs)\n\nTo improve interpretability and stability:\n- Only **electricity meters** are used\n- Modeling is restricted to a **single site**\n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nplt.rcParams[\"figure.figsize\"] = (8, 4)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:31:51.472364Z","iopub.execute_input":"2026-01-20T08:31:51.473112Z","iopub.status.idle":"2026-01-20T08:31:51.477672Z","shell.execute_reply.started":"2026-01-20T08:31:51.473076Z","shell.execute_reply":"2026-01-20T08:31:51.476771Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data Loading\n\nThe notebook assumes the ASHRAE Energy Prediction dataset is attached in Kaggle.\n","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/train.csv\")\nbuilding_df = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/building_metadata.csv\")\nweather_df = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/weather_train.csv\")\n\ntrain_df[\"timestamp\"] = pd.to_datetime(train_df[\"timestamp\"])\nweather_df[\"timestamp\"] = pd.to_datetime(weather_df[\"timestamp\"])\n\nprint(train_df.shape, building_df.shape, weather_df.shape)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:31:51.478803Z","iopub.execute_input":"2026-01-20T08:31:51.479137Z","iopub.status.idle":"2026-01-20T08:32:03.113737Z","shell.execute_reply.started":"2026-01-20T08:31:51.479101Z","shell.execute_reply":"2026-01-20T08:32:03.112644Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Filtering & Scoping\n\nWe restrict the analysis to:\n- Electricity meters only\n- A single site (to avoid climate heterogeneity)\n- Rows with complete critical data\n","metadata":{}},{"cell_type":"code","source":"SITE_ID = 13  # chosen site\n\n# Restrict to electricity meter first\ntrain_site = train_df[train_df[\"meter\"] == 0].copy()\n\n# Merge building metadata (this brings site_id)\ntrain_site = train_site.merge(\n    building_df[[\"building_id\", \"site_id\", \"square_feet\", \"year_built\", \"primary_use\"]],\n    on=\"building_id\",\n    how=\"left\"\n)\n\n# Now filter by site_id (CORRECT PLACE)\ntrain_site = train_site[train_site[\"site_id\"] == SITE_ID].copy()\n\n# Filter weather data for the same site\nweather_site = weather_df[\n    weather_df[\"site_id\"] == SITE_ID\n][[\"site_id\", \"timestamp\", \"air_temperature\", \"dew_temperature\", \"wind_speed\"]].copy()\n\n# Merge weather\ndata_df = train_site.merge(\n    weather_site,\n    on=[\"site_id\", \"timestamp\"],\n    how=\"left\"\n)\n\n# Drop missing critical values\nrequired = [\n    \"meter_reading\",\n    \"air_temperature\",\n    \"dew_temperature\",\n    \"wind_speed\",\n    \"square_feet\"\n]\n\ndata_df = data_df.dropna(subset=required).copy()\n\nprint(\"Final dataset shape:\", data_df.shape)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:03.114859Z","iopub.execute_input":"2026-01-20T08:32:03.115210Z","iopub.status.idle":"2026-01-20T08:32:05.650777Z","shell.execute_reply.started":"2026-01-20T08:32:03.115171Z","shell.execute_reply":"2026-01-20T08:32:05.649920Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Feature Engineering\n\nFeatures:\n- Weather: air temperature, dew temperature, wind speed\n- Time: hour of day, day of week\n- Building scale: log(square_feet)\n\nTarget:\n- log(hourly electricity consumption)\n","metadata":{}},{"cell_type":"code","source":"data_df[\"log_meter_reading\"] = np.log1p(data_df[\"meter_reading\"])\ndata_df[\"log_square_feet\"] = np.log1p(data_df[\"square_feet\"])\n\ndata_df[\"hour\"] = data_df[\"timestamp\"].dt.hour\ndata_df[\"dayofweek\"] = data_df[\"timestamp\"].dt.dayofweek\n\nfeatures = [\n    \"air_temperature\",\n    \"dew_temperature\",\n    \"wind_speed\",\n    \"hour\",\n    \"dayofweek\",\n    \"log_square_feet\"\n]\n\nX = data_df[features].values\ny = data_df[\"log_meter_reading\"].values\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:05.653424Z","iopub.execute_input":"2026-01-20T08:32:05.653731Z","iopub.status.idle":"2026-01-20T08:32:05.815496Z","shell.execute_reply.started":"2026-01-20T08:32:05.653706Z","shell.execute_reply":"2026-01-20T08:32:05.814746Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Time-Aware Train / Validation Split\n\nTraining uses past data.\nValidation uses future data.\n","metadata":{}},{"cell_type":"code","source":"split_time = data_df[\"timestamp\"].quantile(0.8)\n\ntrain_mask = data_df[\"timestamp\"] <= split_time\nval_mask = data_df[\"timestamp\"] > split_time\n\nX_train, X_val = X[train_mask], X[val_mask]\ny_train, y_val = y[train_mask], y[val_mask]\n\nprint(X_train.shape, X_val.shape)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:05.816901Z","iopub.execute_input":"2026-01-20T08:32:05.817234Z","iopub.status.idle":"2026-01-20T08:32:05.891291Z","shell.execute_reply.started":"2026-01-20T08:32:05.817206Z","shell.execute_reply":"2026-01-20T08:32:05.890446Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Feature Normalization\n\nFeatures are normalized using training statistics only.\nThis is critical for LMS stability.\n","metadata":{}},{"cell_type":"code","source":"X_mean = X_train.mean(axis=0)\nX_std = X_train.std(axis=0) + 1e-6\n\nX_train_norm = (X_train - X_mean) / X_std\nX_val_norm = (X_val - X_mean) / X_std\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:05.892331Z","iopub.execute_input":"2026-01-20T08:32:05.892661Z","iopub.status.idle":"2026-01-20T08:32:06.065844Z","shell.execute_reply.started":"2026-01-20T08:32:05.892634Z","shell.execute_reply":"2026-01-20T08:32:06.064987Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Online Learning with LMS (SGD)\n\nThe LMS algorithm updates model parameters one sample at a time\nand converges to a **stochastic equilibrium** on noisy data.\n","metadata":{}},{"cell_type":"code","source":"class LMSRegressor:\n    def __init__(self, n_features, learning_rate=0.01):\n        self.lr = learning_rate\n        self.w = np.zeros(n_features)\n        self.b = 0.0\n        self.loss_history = []\n\n    def predict(self, X):\n        return X @ self.w + self.b\n\n    def update(self, x, y):\n        y_hat = np.dot(self.w, x) + self.b\n        error = y - y_hat\n        self.w += self.lr * error * x\n        self.b += self.lr * error\n        return error\n\n    def fit(self, X, y, n_epochs=5, shuffle=True):\n        n = X.shape[0]\n        for _ in range(n_epochs):\n            idx = np.random.permutation(n) if shuffle else np.arange(n)\n            loss = 0.0\n            for i in idx:\n                e = self.update(X[i], y[i])\n                loss += 0.5 * e**2\n            self.loss_history.append(loss / n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:06.066975Z","iopub.execute_input":"2026-01-20T08:32:06.067279Z","iopub.status.idle":"2026-01-20T08:32:06.075653Z","shell.execute_reply.started":"2026-01-20T08:32:06.067245Z","shell.execute_reply":"2026-01-20T08:32:06.074547Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"lms = LMSRegressor(n_features=X_train_norm.shape[1], learning_rate=0.01)\nlms.fit(X_train_norm, y_train, n_epochs=5)\n\nplt.plot(lms.loss_history, marker=\"o\")\nplt.title(\"LMS Training Loss (Stochastic Equilibrium)\")\nplt.xlabel(\"Epoch\")\nplt.ylabel(\"Loss\")\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:06.076953Z","iopub.execute_input":"2026-01-20T08:32:06.077763Z","iopub.status.idle":"2026-01-20T08:32:32.390750Z","shell.execute_reply.started":"2026-01-20T08:32:06.077725Z","shell.execute_reply":"2026-01-20T08:32:32.389910Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Uncertainty Modeling (Monte Carlo)\n\nWeather features are treated as uncertain.\nTime and building features are deterministic.\n","metadata":{}},{"cell_type":"code","source":"weather_feature_indices = [0, 1, 2]\nfeature_stds = X_train_norm.std(axis=0)\n\nnoise_scales = {i: 0.05 * feature_stds[i] for i in weather_feature_indices}\n\ndef monte_carlo_perturb(x, noise_scales, n_mc=300):\n    n_points, n_features = x.shape\n    X_mc = np.repeat(x[None, :, :], n_mc, axis=0)\n    for idx, scale in noise_scales.items():\n        noise = np.random.normal(0.0, scale, size=(n_mc, n_points))\n        X_mc[:, :, idx] += noise\n    return X_mc\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:32.391889Z","iopub.execute_input":"2026-01-20T08:32:32.392193Z","iopub.status.idle":"2026-01-20T08:32:32.483920Z","shell.execute_reply.started":"2026-01-20T08:32:32.392160Z","shell.execute_reply":"2026-01-20T08:32:32.482852Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Monte Carlo Prediction & Uncertainty Analysis\n","metadata":{}},{"cell_type":"code","source":"N_POINTS = 50\nN_MC = 300\n\nX_batch = X_val_norm[:N_POINTS]\ny_det = lms.predict(X_batch)\n\nX_mc = monte_carlo_perturb(X_batch, noise_scales, n_mc=N_MC)\n\ny_mc = np.array([lms.predict(X_mc[k]) for k in range(N_MC)])\n\ny_mean = y_mc.mean(axis=0)\ny_std = y_mc.std(axis=0)\ny_lo = np.percentile(y_mc, 2.5, axis=0)\ny_hi = np.percentile(y_mc, 97.5, axis=0)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:32.485336Z","iopub.execute_input":"2026-01-20T08:32:32.485738Z","iopub.status.idle":"2026-01-20T08:32:32.497770Z","shell.execute_reply.started":"2026-01-20T08:32:32.485687Z","shell.execute_reply":"2026-01-20T08:32:32.496772Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"idx = 0\n\nplt.hist(y_mc[:, idx], bins=40, alpha=0.7)\nplt.axvline(y_det[idx], color=\"red\", linestyle=\"--\", label=\"Deterministic\")\nplt.axvline(y_mean[idx], color=\"black\", label=\"MC Mean\")\nplt.title(\"Prediction Distribution (Single Sample)\")\nplt.legend()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:32.499230Z","iopub.execute_input":"2026-01-20T08:32:32.499912Z","iopub.status.idle":"2026-01-20T08:32:32.710172Z","shell.execute_reply.started":"2026-01-20T08:32:32.499870Z","shell.execute_reply":"2026-01-20T08:32:32.709062Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.plot(y_det, label=\"Deterministic\", marker=\"o\")\nplt.plot(y_mean, label=\"MC Mean\", marker=\"x\")\nplt.fill_between(range(N_POINTS), y_lo, y_hi, alpha=0.3, label=\"95% CI\")\nplt.title(\"Uncertainty-Aware Predictions\")\nplt.legend()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-20T08:32:32.711347Z","iopub.execute_input":"2026-01-20T08:32:32.711698Z","iopub.status.idle":"2026-01-20T08:32:32.951174Z","shell.execute_reply.started":"2026-01-20T08:32:32.711662Z","shell.execute_reply":"2026-01-20T08:32:32.949970Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Key Takeaways\n\n- LMS converges to a **stochastic equilibrium** on noisy data\n- Point predictions hide risk\n- Monte Carlo simulation reveals **prediction confidence**\n- The same model can produce very different uncertainty levels\n\nThis notebook focuses on **methodology and uncertainty reasoning**, not deployment infrastructure.\n","metadata":{}}]}