{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Bayesian Hyperparameter Optimization\n\nFor competitions as close as those in Kaggle, the amount of lift that one can get from optimizing hyperparameters often becomes the decisive factor. It seems to be even more the case in this competition where we have to predict future prices of something as volatile as Crypto. \n\n# Data Loading and All That Stuff\n\nCheck out this amazing notebook! \n\n[https://www.kaggle.com/craniket/gresearch-submitting-lagged-features-via-api?scriptVersionId=85838940](http://)","metadata":{}},{"cell_type":"code","source":"\nimport os\nimport random\nimport pandas as pd\nimport numpy as np\nimport lightgbm as lgb\nfrom lightgbm import LGBMRegressor\nimport gresearch_crypto\nimport time\nimport datetime\nfrom hyperopt import tpe\nfrom hyperopt import hp\nfrom hyperopt import Trials, STATUS_OK,fmin\nfrom hyperopt.pyll.base import scope\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.metrics import make_scorer\nTRAIN_CSV = '/kaggle/input/g-research-crypto-forecasting/train.csv'\nASSET_DETAILS_CSV = '/kaggle/input/g-research-crypto-forecasting/asset_details.csv'","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:38:39.545081Z","iopub.execute_input":"2022-01-20T06:38:39.545461Z","iopub.status.idle":"2022-01-20T06:38:42.606073Z","shell.execute_reply.started":"2022-01-20T06:38:39.545372Z","shell.execute_reply":"2022-01-20T06:38:42.605347Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train = pd.read_csv(TRAIN_CSV)\ndf_train.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:38:42.607745Z","iopub.execute_input":"2022-01-20T06:38:42.607975Z","iopub.status.idle":"2022-01-20T06:39:36.312982Z","shell.execute_reply.started":"2022-01-20T06:38:42.607943Z","shell.execute_reply":"2022-01-20T06:39:36.312316Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_asset_details = pd.read_csv(ASSET_DETAILS_CSV).sort_values(\"Asset_ID\")\ndf_asset_details","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:36.314197Z","iopub.execute_input":"2022-01-20T06:39:36.314517Z","iopub.status.idle":"2022-01-20T06:39:36.335245Z","shell.execute_reply.started":"2022-01-20T06:39:36.314479Z","shell.execute_reply":"2022-01-20T06:39:36.334413Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_features(df, \n                 asset_id, \n                 train=True):\n    '''\n    This function takes a dataframe with all asset data and return the lagged features for a single asset.\n    \n    df - Full dataframe with all assets included\n    asset_id - integer from 0-13 inclusive to represent a cryptocurrency asset\n    train - True - you are training your model\n          - False - you are submitting your model via api\n    '''\n    \n    df = df[df['Asset_ID']==asset_id]\n    df = df.sort_values('timestamp')\n    if train == True:\n        df_feat = df.copy()\n        # define a train_flg column to split your data into train and validation\n        totimestamp = lambda s: np.int32(time.mktime(datetime.datetime.strptime(s, \"%d/%m/%Y\").timetuple()))\n        valid_window = [totimestamp(\"12/03/2021\")]\n        df_feat['train_flg'] = np.where(df_feat['timestamp']>=valid_window[0], 0,1)\n        df_feat = df_feat[['timestamp','Asset_ID','Close','Target','train_flg']].copy()\n    else:\n        df = df.sort_values('row_id')\n        df_feat = df[['Asset_ID','Close','row_id']].copy()\n    \n    # Create your features here, they can be lagged or not\n    df_feat['sma15'] = df_feat['Close'].rolling(15).mean()/df_feat['Close'] -1\n    df_feat['sma60'] = df_feat['Close'].rolling(60).mean()/df_feat['Close'] -1\n    df_feat['sma240'] = df_feat['Close'].rolling(240).mean()/df_feat['Close'] -1\n    \n    df_feat['return15'] = df_feat['Close']/df_feat['Close'].shift(15) -1\n    df_feat['return60'] = df_feat['Close']/df_feat['Close'].shift(60) -1\n    df_feat['return240'] = df_feat['Close']/df_feat['Close'].shift(240) -1\n    df_feat = df_feat.fillna(0)\n    \n    return df_feat","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:36.337567Z","iopub.execute_input":"2022-01-20T06:39:36.3379Z","iopub.status.idle":"2022-01-20T06:39:36.349295Z","shell.execute_reply.started":"2022-01-20T06:39:36.337861Z","shell.execute_reply":"2022-01-20T06:39:36.348094Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create your feature dataframe for each asset and concatenate\nfeature_df = pd.DataFrame()\nfor i in range(14):\n    feature_df = pd.concat([feature_df,get_features(df_train,i,train=True)])","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:36.35103Z","iopub.execute_input":"2022-01-20T06:39:36.351305Z","iopub.status.idle":"2022-01-20T06:39:53.599465Z","shell.execute_reply.started":"2022-01-20T06:39:36.351263Z","shell.execute_reply":"2022-01-20T06:39:53.598737Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# assign weight column feature dataframe\nfeature_df = pd.merge(feature_df, df_asset_details[['Asset_ID','Weight']], how='left', on=['Asset_ID'])","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:53.601343Z","iopub.execute_input":"2022-01-20T06:39:53.602118Z","iopub.status.idle":"2022-01-20T06:39:56.775608Z","shell.execute_reply.started":"2022-01-20T06:39:53.60208Z","shell.execute_reply":"2022-01-20T06:39:56.774872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# define features for LGBM\nfeatures = ['Asset_ID','sma15','sma60','sma240','return15','return60','return240','Weight']\ncategoricals = ['Asset_ID']","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:56.776719Z","iopub.execute_input":"2022-01-20T06:39:56.776958Z","iopub.status.idle":"2022-01-20T06:39:56.781214Z","shell.execute_reply.started":"2022-01-20T06:39:56.776925Z","shell.execute_reply":"2022-01-20T06:39:56.78053Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del df_train","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:56.782437Z","iopub.execute_input":"2022-01-20T06:39:56.782851Z","iopub.status.idle":"2022-01-20T06:39:56.795991Z","shell.execute_reply.started":"2022-01-20T06:39:56.782814Z","shell.execute_reply":"2022-01-20T06:39:56.795151Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# define the evaluation metric\ndef weighted_correlation(a, train_data):\n    \n    weights = train_data.Weight.values\n    b = train_data.target\n    \n    \n    w = np.ravel(weights)\n    a = np.ravel(a)\n    b = np.ravel(b)\n\n    sum_w = np.sum(w)\n    mean_a = np.sum(a * w) / sum_w\n    mean_b = np.sum(b * w) / sum_w\n    var_a = np.sum(w * np.square(a - mean_a)) / sum_w\n    var_b = np.sum(w * np.square(b - mean_b)) / sum_w\n\n    cov = np.sum((a * b * w)) / np.sum(w) - mean_a * mean_b\n    corr = cov / np.sqrt(var_a * var_b)\n\n    return corr","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:56.79721Z","iopub.execute_input":"2022-01-20T06:39:56.797498Z","iopub.status.idle":"2022-01-20T06:39:56.806824Z","shell.execute_reply.started":"2022-01-20T06:39:56.797446Z","shell.execute_reply":"2022-01-20T06:39:56.806055Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Tuning\n\nHere is the part where we use Bayesian hyperparameter tuning. Simply put, bayesian hyperparameter tuning is an \"intelligent\" mechanism wherein each model evaluation carries information from the previous iterations. Compare this to gridsearch or random search, which is an [embarassingly parallel](https://en.wikipedia.org/wiki/Embarrassingly_parallel) problem. Both of these methods involve evaluating the objective function at a bunch of points on the hyperparameter space separately from each other. However, note that we used the term \"iterations\" for Bayesian hyperparameter tuning.\n\n# How does it work?\n# What does hyperparameter optimization entail?\n\nWhat do we REALLY want to find when we optimize hyperparameters? We want to find the set of values for which the loss function or objective function is minimized. The objective function depends on \n* The observations\n* The hyperparameters\n\n\nWhat happens if we fix the hyperparameters? The loss function becomes analogous to probability distributions - it takes up values in ranges depending on the distribution of the predictor variables.\n\nThis can be denoted by ***P(loss|hyperparameters)***.\n\nOnce we have this, we can estimate the loss for this set of hyperparameters with something like the  median of the assumed distribution.\n\nHowever, we do not know this distribution for \"each\" value of the hyperparameter. The idea is to figure out this distribution at a sufficient number of points. That is, sufficient to get a decent idea of the joint distribution. \n\n# How many points is sufficient?\n\nLet's try to understand this visually. I am shamelessly stealing these pictures from [this amazing lecture](https://www.cs.toronto.edu/~rgrosse/courses/csc321_2017/slides/lec21.pdf).\n\nLet's say we want to approximate a function, and we have evaluated its value at three points.\n\n\n![](https://i.imgur.com/ujmYbRB.png)\n\nOne non-linear estimate of the function could be\n\n\n![](https://i.imgur.com/hKeBPjt.png)\n\nNote the labelled 80, 90, 95th percentiles. We are extremely certain of the function's value at points we have evaluated it. We grow less and less certain as we move away.\n\nWe want to evaluate the function at enough points so that we are \"pretty\" certain of the function. And here is where the Bayesian optimization process comes in - it tells us where next to evaluate the function so as to gain most certainty. \n\nLet's look at the components here.\n\n# Components \n\nWe obviously have the objective function to think about. However, in practice, it is too complicated to estimate/approximate. People tend to use different kinds of gaussian kernels as proxies, or **surrogate functions**.\n\nAnd then there's the function which is going to tell us which point to look at next. This is called the **acquisition function**. In practice, people tend to use the expected  improvement, E[max(γ − f (θ), 0)], where γ is the previous minima found. The θ at which this is maximum is chosen for the next evaluation. This is how it could look:\n\n\n![](https://i.imgur.com/OAXvXTj.png)\n\nIn this notebook, we will use the TPE algorithm, which uses gaussian kernels and expected improvement. It uses Bayes Theorem to calculate the expected improvement. The expectation is calculated over the conditional surrogate function, or over *P(loss|hyperparameters)*. TPE expresses *P(loss|hyperparameters)* as a function of *P(hyperparameters|loss)*. Now *P(hyperparameters|loss)* is defined as \n\n\n![](https://i.imgur.com/LcHJN92.png)\n\n\nwhere y* is a previously defined threshold. This creates two distributions, one for scores lesser than the threshold and the other for scores higher than the threshold. The idea is to make it more likely to draw hyperparameters from regions resulting in lower losses. And indeed, when the expectation is calculated, it is found to be proportional to the ratio *l(x)/g(x)*!\n\nMore details can be found in this brilliant article https://towardsdatascience.com/a-conceptual-explanation-of-bayesian-model-based-hyperparameter-optimization-for-machine-learning-b8172278050f\n\nAnd finally, there is, of course, the hyperparameter space over which we want to find the minima.\n\nTo summarize, the components are:\n\n* Surrogate function\n* Acquisition function\n* Hyperparameter space\n\nLet's execute all of this. We will use the hyperopt package in python.\n","metadata":{}},{"cell_type":"code","source":"def hyperopt(param_space, X_train, y_train, X_test, y_test, num_eval):\n\n    start = time.time()\n    \n    def objective_function(params):\n        clf = LGBMRegressor(**params).fit(X_train,y_train)\n        y_pred = clf.predict(X_test)\n        X_test.target = y_test\n        score = weighted_correlation(y_pred,X_test)\n        return {'loss': -score, 'status': STATUS_OK}\n\n    trials = Trials()\n    best_param = fmin(objective_function, \n                      param_space, \n                      algo=tpe.suggest, \n                      max_evals=num_eval, \n                      trials=trials)\n    loss = [x['result']['loss'] for x in trials.trials]\n    \n    best_param_values = [x for x in best_param.values()]\n\n    print(\"\")\n    print(\"##### Results\")\n    print(\"Best parameters: \", best_param)\n    print(\"Best score: \", -min(loss))\n    print(\"Time elapsed: \", time.time() - start)\n    print(\"Parameter combinations evaluated: \", num_eval)\n    \n    return trials","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:56.809184Z","iopub.execute_input":"2022-01-20T06:39:56.809488Z","iopub.status.idle":"2022-01-20T06:39:56.820191Z","shell.execute_reply.started":"2022-01-20T06:39:56.809454Z","shell.execute_reply":"2022-01-20T06:39:56.819471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# define train and validation weights and datasets\n#weights_train = feature_df.query('train_flg == 1')[['Weight']]\n#weights_test = feature_df.query('train_flg == 0')[['Weight']]\n\n#train_dataset = lgb.Dataset(feature_df.query('train_flg == 1')[features], \n#                            feature_df.query('train_flg == 1')['Target'].values, \n#                            feature_name = features, \n#                            categorical_feature= categoricals)\n#val_dataset = lgb.Dataset(feature_df.query('train_flg == 0')[features], \n#                          feature_df.query('train_flg == 0')['Target'].values, \n#                          feature_name = features, \n#                          categorical_feature= categoricals)\n\n#train_dataset.add_w = weights_train\n#val_dataset.add_w = weights_test\n\nevals_result = {}\n\nspace = { 'learning_rate': hp.loguniform('learning_rate', np.log(0.01), np.log(0.2)),\n    'max_depth': scope.int(hp.quniform('max_depth', 4, 10, 1)),\n    'n_estimators': scope.int(hp.quniform('n_estimators', 200, 2000, 50)),\n#    'early_stopping_rounds': scope.int(hp.quniform('early_stopping_rounds', 20, 500, 5)),\n    'objective': 'regression',\n    'metric': 'None',\n    'boosting_type': hp.choice('boosting_type',['goss','gbdt']),\n    'verbose': -1,\n    'seed': 46\n}\n\nresults_hyperopt = hyperopt(space, feature_df.query('train_flg == 1')[features], feature_df.query('train_flg == 1')['Target'].values, feature_df.query('train_flg == 0')[features], feature_df.query('train_flg == 0')['Target'].values, 25)","metadata":{"execution":{"iopub.status.busy":"2022-01-20T06:39:56.82317Z","iopub.execute_input":"2022-01-20T06:39:56.823382Z","iopub.status.idle":"2022-01-20T14:08:48.17354Z","shell.execute_reply.started":"2022-01-20T06:39:56.823358Z","shell.execute_reply":"2022-01-20T14:08:48.172717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### Results\n\nBest parameters:  {'boosting_type': 1, 'learning_rate': 0.010257600094741514, 'max_depth': 4.0, 'n_estimators': 1550.0}\nBest score:  0.03316784490207032\nTime elapsed:  26925.81499361992\nParameter combinations evaluated:  25","metadata":{"_kg_hide-input":true}},{"cell_type":"markdown","source":"Thanks for checking this out! Let me know what you think! ","metadata":{}}]}