{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Target Encoding \n\n## Target Encoding in Practice\n\nWe encode categorical variables using the mean of the target variable within a category/group; for example, \n```\ntarget_col = \"meter_reading\"\ngroup_cols = [\"building_id\", \"meter\"]\ntarget_encoding = df.groupby(group_cols)[target_col].transform(\"mean\")\n```\n\n## Target Encoding in Theory\nWe treat the target variable $y$ within each category as a random variable.  The distribution of the random variable is problem-dependent. In this competition, we assume that the meter readings are log-normal; in other words, `log1p(meter_reading)` is normaly distributed.  Target encoding can be viewed as estimating the parameters of the normal distribution.  \n\n**Frequentist Approach**: The frequentist approach to parameter estimation is maximum likelihood estimation (mle), which finds the parameters $\\theta = (\\mu,\\sigma)$ that maximize the likelihood function $f(y \\,\\vert\\, \\theta)$.  For a normal distribution, this results in the following maximum likelihood estimates:\n\nthe mean\n$$\n\\mu_{mle} = \\frac{1}{n}\\sum_{i=1}^n y[i]\n$$\n\nthe variance\n$$\n\\sigma^2_{mle} = \\frac{1}{n}\\sum_{i=1}^n (y[i] - \\mu_{mle})^2.\n$$\n\n**Bayesian Approach**: The Bayesian approach assumes a prior distribution $f(\\theta)$ over the parameters of the distribution and uses data and bayes rule to calculate an updated (the  posterior) distribution \n\n$$\nf(\\theta \\,\\vert\\, y) = \\frac{f(y \\,\\vert\\, \\theta) \\, f(\\theta)}{f(y)}.\n$$\n\nFor simplicity, we assume the precision of the likelihood is known and equal to the one over the maximum likelihood estimate $\\tau = 1 / \\sigma^2_{mle}$. This leaves the mean $\\mu$ as the only parameter to be estimated. We take the mean of the prior distribution $\\mu^2_{prior}$ to be the global average meter reading. We treat the precision of the prior distribution $\\tau_{prior} = 1/\\sigma^2_{prior}$ as a regularizing hyperparameter that can be tuned using cross-validation.  Larger values of $\\tau_{prior}$ will make the posterior estimates closer to their prior counterparts; in other words, increasing $\\tau_{prior}$ will increase the amount of regularization. The parameters of the posterior distribution are given by the following equations:\n\nthe mean\n$$\n\\mu_{post} = \\frac {\n    \\tau_{prior} \\mu _{prior} +\n    n \\tau \\mu_{mle}\n}{\n    \\tau_{prior} +  \n    n \\tau\n}\n$$\n\nthe precision\n$$\n\\tau_{post} = \\tau_{prior}+n\\tau.\n$$\n\nTo gain some intution about the posterior mean, note that it is a weighted average of the prior mean and the maximum likelihood estiamte. To see this, let \n$$\n\\alpha = \\frac{\\tau_{prior}}{\\tau_{prior} +  n \\tau} \n$$ \nthus\n$$\n\\mu_{post} = \\alpha \\mu_{prior}  + (1 - \\alpha) \\mu_{mle}.\n$$\n\nThe most natural way to encode the category is using the posterior mean, but it also makes sense to encode any other statistic of the posterior\ndistribution. For example, the posterior precision would tell us how certain we are about the estimate for the category's target variable mean. \n\n\n**References**: \n* [Wikipedia page on Conjugate Priors](https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution)\n* [MIT 18.05 Slides on Conjugate Priors](https://math.mit.edu/~dav/05.dir/class15-slides-all.pdf)\n\n\n## Hierachical Bayesian Target Encoding\n\nIn the previous section we used the global average meter reading as the prior mean. However, we could take advantage of the hierarchy in the data to define a more informative prior. For example, we could use the mean target encoded value for `meter` as the prior mean for the `building_id`, `meter` grouping. "},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import gc\nimport numpy as np\nimport pandas as pd \nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.metrics import mean_squared_error","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"PRIOR_PRECISION = 10","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class GaussianTargetEncoder():\n        \n    def __init__(self, group_cols, target_col=\"target\", prior_cols=None):\n        self.group_cols = group_cols\n        self.target_col = target_col\n        self.prior_cols = prior_cols\n\n    def _get_prior(self, df):\n        if self.prior_cols is None:\n            prior = np.full(len(df), df[self.target_col].mean())\n        else:\n            prior = df[self.prior_cols].mean(1)\n        return prior\n                    \n    def fit(self, df):\n        self.stats = df.assign(mu_prior=self._get_prior(df), y=df[self.target_col])\n        self.stats = self.stats.groupby(self.group_cols).agg(\n            n        = (\"y\", \"count\"),\n            mu_mle   = (\"y\", np.mean),\n            sig2_mle = (\"y\", np.var),\n            mu_prior = (\"mu_prior\", np.mean),\n        )        \n    \n    def transform(self, df, prior_precision=1000, stat_type=\"mean\"):\n        \n        precision = prior_precision + self.stats.n/self.stats.sig2_mle\n        \n        if stat_type == \"mean\":\n            numer = prior_precision*self.stats.mu_prior\\\n                    + self.stats.n/self.stats.sig2_mle*self.stats.mu_mle\n            denom = precision\n        elif stat_type == \"var\":\n            numer = 1.0\n            denom = precision\n        elif stat_type == \"precision\":\n            numer = precision\n            denom = 1.0\n        else: \n            raise ValueError(f\"stat_type={stat_type} not recognized.\")\n        \n        mapper = dict(zip(self.stats.index, numer / denom))\n        if isinstance(self.group_cols, str):\n            keys = df[self.group_cols].values.tolist()\n        elif len(self.group_cols) == 1:\n            keys = df[self.group_cols[0]].values.tolist()\n        else:\n            keys = zip(*[df[x] for x in self.group_cols])\n        \n        values = np.array([mapper.get(k) for k in keys]).astype(float)\n        \n        prior = self._get_prior(df)\n        values[~np.isfinite(values)] = prior[~np.isfinite(values)]\n        \n        return values\n    \n    def fit_transform(self, df, *args, **kwargs):\n        self.fit(df)\n        return self.transform(df, *args, **kwargs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def rmsle(x,y):\n    x = np.log1p(x)\n    y = np.log1p(y)\n    return np.sqrt(mean_squared_error(x, y))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load data\ntrain = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/train.csv\")\ntest  = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/test.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# sample for kernel\ntrain = train.sample(int(2.5e5)).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create target\ntrain[\"target\"] = np.log1p(train.meter_reading)\ntest[\"target\"] = train.target.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create time features\ndef add_time_features(df):\n    df.timestamp = pd.to_datetime(df.timestamp)    \n    df[\"hour\"]    = df.timestamp.dt.hour\n    df[\"weekday\"] = df.timestamp.dt.weekday\n    df[\"month\"]   = df.timestamp.dt.month\n\nadd_time_features(train)\nadd_time_features(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# define groupings and corresponding priors\ngroups_and_priors = {\n    \n    # singe encodings\n    (\"hour\",):        None,\n    (\"weekday\",):     None,\n    (\"month\",):       None,\n    (\"building_id\",): None,\n    (\"meter\",):       None,\n    \n    # second-order interactions\n    (\"meter\", \"hour\"):        [\"gte_meter\", \"gte_hour\"],\n    (\"meter\", \"weekday\"):     [\"gte_meter\", \"gte_weekday\"],\n    (\"meter\", \"month\"):       [\"gte_meter\", \"gte_month\"],\n    (\"meter\", \"building_id\"): [\"gte_meter\", \"gte_building_id\"],\n        \n    # higher-order interactions\n    (\"meter\", \"building_id\", \"hour\"):    [\"gte_meter_building_id\", \"gte_meter_hour\"],\n    (\"meter\", \"building_id\", \"weekday\"): [\"gte_meter_building_id\", \"gte_meter_weekday\"],\n    (\"meter\", \"building_id\", \"month\"):   [\"gte_meter_building_id\", \"gte_meter_month\"],\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"features = []\nfor group_cols, prior_cols in groups_and_priors.items():\n    features.append(f\"gte_{'_'.join(group_cols)}\")\n    gte = GaussianTargetEncoder(list(group_cols), \"target\", prior_cols)    \n    train[features[-1]] = gte.fit_transform(train, PRIOR_PRECISION)\n    test[features[-1]]  = gte.transform(test,  PRIOR_PRECISION)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# clean up\ndrop_cols = [\"hour\", \"weekday\", \"month\", \"building_id\"]\ntrain.drop(drop_cols, 1, inplace=True)\ntest.drop(drop_cols, 1, inplace=True)\ndel  gte\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train[features + [\"target\"]].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test[features].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_preds = np.zeros(len(train))\ntest_preds = np.zeros(len(test))\n\nfor m in range(4):\n    \n    print(f\"Meter {m}\", end=\"\") \n    \n    # instantiate model\n    model = RidgeCV(\n        alphas=np.logspace(-10, 1, 25), \n        normalize=True,\n    )    \n    \n    # fit model\n    model.fit(\n        X=train.loc[train.meter==m, features].values, \n        y=train.loc[train.meter==m, \"target\"].values\n    )\n\n    # make predictions \n    train_preds[train.meter==m] = model.predict(train.loc[train.meter==m, features].values)\n    test_preds[test.meter==m]   = model.predict(test.loc[test.meter==m, features].values)\n    \n    # transform predictions\n    train_preds[train_preds < 0] = 0\n    train_preds[train.meter==m] = np.expm1(train_preds[train.meter==m])\n    \n    test_preds[test_preds < 0] = 0 \n    test_preds[test.meter==m] = np.expm1(test_preds[test.meter==m])\n    \n    # evaluate model\n    meter_rmsle = rmsle(\n        train_preds[train.meter==m],\n        train.loc[train.meter==m, \"meter_reading\"].values\n    )\n    \n    print(f\", rmsle={meter_rmsle:0.5f}\")\n\nprint(f\"Overall rmsle={rmsle(train_preds, train.meter_reading.values):0.5f}\")\ndel train, train_preds, test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create submission\nsubm  = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/sample_submission.csv\")\nsubm[\"meter_reading\"] = test_preds\nsubm.to_csv(f\"submission.csv\", index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.5"}},"nbformat":4,"nbformat_minor":1}