{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.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":21651,"databundleVersionId":1595136,"sourceType":"competition"},{"sourceId":6965130,"sourceType":"datasetVersion","datasetId":4001431},{"sourceId":6966996,"sourceType":"datasetVersion","datasetId":4002702}],"dockerImageVersionId":30579,"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","execution":{"iopub.status.busy":"2024-04-28T11:36:35.467744Z","iopub.execute_input":"2024-04-28T11:36:35.468149Z","iopub.status.idle":"2024-04-28T11:36:35.961435Z","shell.execute_reply.started":"2024-04-28T11:36:35.468111Z","shell.execute_reply":"2024-04-28T11:36:35.960182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"questions = pd.read_csv('../input/riiid-test-answer-prediction/questions.csv',                         \n                        usecols=[0, 1, 2, 3, 4],\n                           dtype={'question_id': 'int16',\n                                  'part': 'int8',\n                                  'bundle_id': 'int8',\n                                  'correct_answer': 'int8',\n                                  'tags': 'object'}\n                          )\n\n\ndf_train = pd.read_csv('../input/riiid-test-answer-prediction/train.csv',\n                   usecols=[0, 1, 2, 3, 4, 5, 7, 8, 9],\n                   dtype={\n                          'row_id': 'int32',\n                       'timestamp': 'int64',\n                          'user_id': 'int32',\n                          'content_id': 'int16',\n                          'content_type_id': 'int8',\n                          'task_container_id': 'int16',\n                          'user_answer': 'int8',\n                          'answered_correctly':'int8',\n                          'prior_question_elapsed_time': 'float32',\n                          'prior_question_had_explanation': 'boolean'\n                          },\n                   )\n\nprint(df_train.shape)\n\ndf_train = df_train[df_train.answered_correctly!= -1 ]\ndf_train = df_train[df_train['content_type_id'] == 0]\n\ndf_train.drop(['content_type_id'], axis=1, inplace=True)\n\nprint(df_train.shape)\n\ndf_train.dropna(inplace=True)\n\ndf_train = df_train.drop_duplicates(subset=['user_id', 'content_id'], keep='first')\n\nprint(df_train.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T11:36:35.964060Z","iopub.execute_input":"2024-04-28T11:36:35.965075Z","iopub.status.idle":"2024-04-28T11:41:51.368665Z","shell.execute_reply.started":"2024-04-28T11:36:35.965027Z","shell.execute_reply":"2024-04-28T11:41:51.366476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sort rows by timestamp in each item\ndf_train['prior_question_elapsed_time_shift'] = df_train.sort_values(by = 'timestamp').groupby('user_id')['prior_question_elapsed_time'].shift(-1)\ndf_train.dropna(inplace=True)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T11:41:51.372322Z","iopub.execute_input":"2024-04-28T11:41:51.372854Z","iopub.status.idle":"2024-04-28T11:43:12.627280Z","shell.execute_reply.started":"2024-04-28T11:41:51.372811Z","shell.execute_reply":"2024-04-28T11:43:12.626007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#remove outliers for each question\nimport tqdm\n\ntime_upper_bound = df_train.groupby(['content_id'])['prior_question_elapsed_time_shift'].quantile(0.99)\ntime_lower_bound = df_train.groupby(['content_id'])['prior_question_elapsed_time_shift'].quantile(0.01)\n\ninclude = []\nfor row in tqdm.tqdm(df_train.itertuples(),total=df_train.shape[0]):\n\n    if row.timestamp <= time_upper_bound[row.content_id] and row.timestamp >= time_lower_bound[row.content_id]:\n        include.append(True)\n    else:\n        include.append(False)\n\ndf_train['include'] = include\ndf_train = df_train[df_train['include']==True]\n\ndf_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-28T11:43:12.630448Z","iopub.execute_input":"2024-04-28T11:43:12.630841Z","iopub.status.idle":"2024-04-28T11:59:44.351359Z","shell.execute_reply.started":"2024-04-28T11:43:12.630808Z","shell.execute_reply":"2024-04-28T11:59:44.349781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## keep sample size >20 \ntmp = df_train['content_id'].value_counts().reset_index()\nquestion_id_gt_20 = tmp[tmp['count'] > 20]['content_id'].tolist()\ndf_train = df_train[df_train['content_id'].isin(question_id_gt_20)]\ndf_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-28T11:59:44.352997Z","iopub.execute_input":"2024-04-28T11:59:44.353375Z","iopub.status.idle":"2024-04-28T11:59:44.443766Z","shell.execute_reply.started":"2024-04-28T11:59:44.353343Z","shell.execute_reply":"2024-04-28T11:59:44.441933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_train['content_id'].value_counts(), df_train['user_id'].value_counts())","metadata":{"execution":{"iopub.status.busy":"2024-04-28T11:59:44.445436Z","iopub.execute_input":"2024-04-28T11:59:44.445819Z","iopub.status.idle":"2024-04-28T11:59:44.529097Z","shell.execute_reply.started":"2024-04-28T11:59:44.445790Z","shell.execute_reply":"2024-04-28T11:59:44.528212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prepare data for LMM model\nPopn1 = df_train.groupby('content_id', sort=False).size().reset_index()\nPopn1.to_csv('/kaggle/working/Popn1.csv')\n\nitem_diff = df_train.groupby('content_id', sort=False, as_index=False)['answered_correctly'].mean().reset_index()\ndf_train = df_train.merge(item_diff, how='inner', left_on='content_id', right_on='content_id')\ndata_pop = df_train.rename(columns = {'answered_correctly_y': 'item_diff'})\n\nXmean_df = data_pop.groupby('content_id', sort=False)[['item_diff']].mean().reset_index()\nXmean_df.to_csv('/kaggle/working/Xmean1.csv')","metadata":{"execution":{"iopub.status.busy":"2024-04-28T11:59:44.530542Z","iopub.execute_input":"2024-04-28T11:59:44.531964Z","iopub.status.idle":"2024-04-28T11:59:44.902078Z","shell.execute_reply.started":"2024-04-28T11:59:44.531870Z","shell.execute_reply":"2024-04-28T11:59:44.900593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# split train, val and test data\n\nfraction = 'first_half'\ntrain_ = df_train.groupby('index').apply(lambda x: x.iloc[:x.content_id.size//2]) #the first half of rows within each item\nval_ = df_train.groupby('index').apply(lambda x: x.iloc[x.content_id.size//2:(x.content_id.size//2 + x.content_id.size//4)]) #the first half of rows within each item\ntest_ = df_train.groupby('index').apply(lambda x: x.iloc[(x.content_id.size//2 + x.content_id.size//4):]) #the first half of rows within each item\n\n\n# train_.to_csv('/kaggle/working/train_{}.csv'.format(fraction))\n# val_.to_csv('/kaggle/working/val_{}.csv'.format(fraction))\n# test_.to_csv('/kaggle/working/test_{}.csv'.format(fraction))\n# train_.groupby('content_id')['prior_question_elapsed_time_shift'].quantile(q=0.95).reset_index().to_csv('/kaggle/working/time_95_{}.csv'.format(fraction))","metadata":{"execution":{"iopub.status.busy":"2024-04-28T11:59:44.903949Z","iopub.execute_input":"2024-04-28T11:59:44.904407Z","iopub.status.idle":"2024-04-28T12:00:01.672174Z","shell.execute_reply.started":"2024-04-28T11:59:44.904366Z","shell.execute_reply":"2024-04-28T12:00:01.670094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_['content_id'].value_counts(), train_['user_id'].value_counts().size)\nprint(val_['content_id'].value_counts(), val_['user_id'].value_counts().size)\nprint(test_['content_id'].value_counts(), test_['user_id'].value_counts().size)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:01.674285Z","iopub.execute_input":"2024-04-28T12:00:01.674805Z","iopub.status.idle":"2024-04-28T12:00:01.779769Z","shell.execute_reply.started":"2024-04-28T12:00:01.674770Z","shell.execute_reply":"2024-04-28T12:00:01.778665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### regression","metadata":{}},{"cell_type":"code","source":"train_full = train_.merge(questions, how='inner', left_on='content_id', right_on='question_id')\nval_full = val_.merge(questions, how='inner', left_on='content_id', right_on='question_id')\ntest_full = test_.merge(questions, how='inner', left_on='content_id', right_on='question_id')","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:01.784955Z","iopub.execute_input":"2024-04-28T12:00:01.786517Z","iopub.status.idle":"2024-04-28T12:00:01.894935Z","shell.execute_reply.started":"2024-04-28T12:00:01.786356Z","shell.execute_reply":"2024-04-28T12:00:01.893712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## find the 95th RT in each item and use it as dependent variable in regression, merge with question feature\n\ntrain_reg = train_.groupby('content_id')['prior_question_elapsed_time_shift'].quantile(q=0.95).reset_index()\ntest_reg_n = train_.groupby('content_id').size().reset_index() #add n\ntrain_reg = train_reg.merge(test_reg_n, how='inner', on='content_id')\n\ntrain_reg = train_reg.merge(questions, how='inner', left_on='content_id', right_on='question_id')\ntrain_reg.rename(columns = {'prior_question_elapsed_time_shift': 'time_q95', 0: 'n'}, inplace = True)\n\nval_reg = val_.groupby('content_id')['prior_question_elapsed_time_shift'].quantile(q=0.95).reset_index()\nval_reg = val_reg.merge(train_reg[['content_id', 'n']], how='inner', on='content_id' )\nval_reg = val_reg.merge(questions, how='inner', left_on='content_id', right_on='question_id')\nval_reg.rename(columns = {'prior_question_elapsed_time_shift': 'time_q95'}, inplace = True)\n\ntest_reg = test_.groupby('content_id')['prior_question_elapsed_time_shift'].quantile(q=0.95).reset_index()\ntest_reg = test_reg.merge(train_reg[['content_id', 'n']], how='inner', on='content_id' )\ntest_reg = test_reg.merge(questions, how='inner', left_on='content_id', right_on='question_id')\ntest_reg.rename(columns = {'prior_question_elapsed_time_shift': 'time_q95'}, inplace = True)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:01.896602Z","iopub.execute_input":"2024-04-28T12:00:01.897041Z","iopub.status.idle":"2024-04-28T12:00:02.078965Z","shell.execute_reply.started":"2024-04-28T12:00:01.897003Z","shell.execute_reply":"2024-04-28T12:00:02.077193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_reg.shape, test_reg.shape, val_reg.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:02.081512Z","iopub.execute_input":"2024-04-28T12:00:02.082089Z","iopub.status.idle":"2024-04-28T12:00:02.091421Z","shell.execute_reply.started":"2024-04-28T12:00:02.082042Z","shell.execute_reply":"2024-04-28T12:00:02.090196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_reg.to_csv('/kaggle/working/train_reg_{}.csv'.format(fraction))\ntest_reg.to_csv('/kaggle/working/test_reg_{}.csv'.format(fraction))\nval_reg.to_csv('/kaggle/working/val_reg_{}.csv'.format(fraction))","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:02.093010Z","iopub.execute_input":"2024-04-28T12:00:02.094004Z","iopub.status.idle":"2024-04-28T12:00:02.210624Z","shell.execute_reply.started":"2024-04-28T12:00:02.093963Z","shell.execute_reply":"2024-04-28T12:00:02.209055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# test_reg = pd.read_csv('/kaggle/input/data-regg/test_reg_{}.csv'.format(fraction))\n\n# #test_reg = test_reg.iloc[int(test_reg.shape[0]//2):,]\n# val_reg = pd.read_csv('/kaggle/input/data-regg/val_reg_{}.csv'.format(fraction))\n# #val_reg = test_reg.iloc[:int(test_reg.shape[0]//2),]\n\n# train_reg = pd.read_csv('/kaggle/input/data-regg/train_reg_{}.csv'.format(fraction))\n# Xmean1 = pd.read_csv('/kaggle/input/data-regg/Xmean1.csv')\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:02.213982Z","iopub.execute_input":"2024-04-28T12:00:02.215781Z","iopub.status.idle":"2024-04-28T12:00:02.222821Z","shell.execute_reply.started":"2024-04-28T12:00:02.215703Z","shell.execute_reply":"2024-04-28T12:00:02.221120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# include item difficulty variable\nXmean1 = pd.read_csv('/kaggle/working/Xmean1.csv') #Xmean_df above\n\ntrain_reg = train_reg.merge(Xmean1, how='inner', left_on='content_id', right_on='content_id')\nval_reg = val_reg.merge(Xmean1, how='inner', left_on='content_id', right_on='content_id')\ntest_reg = test_reg.merge(Xmean1, how='inner', left_on='content_id', right_on='content_id')\n\ntrain_reg.shape, test_reg.shape, val_reg.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:02.225070Z","iopub.execute_input":"2024-04-28T12:00:02.226025Z","iopub.status.idle":"2024-04-28T12:00:02.260128Z","shell.execute_reply.started":"2024-04-28T12:00:02.225978Z","shell.execute_reply":"2024-04-28T12:00:02.259007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## creating tag feature for regression\nall_tags = []\nfor t in train_reg['tags'].tolist():\n    all_tags.extend(t.split())\n    \nall_tags = np.unique(all_tags)\ndef create_tag_df(df):\n    \n    tags_array = np.array([[0] * len(all_tags)] * df.shape[0])\n    for r, tags in enumerate(df['tags'].tolist()):\n        tags_array[r, [i for i, e in enumerate(all_tags) if e in tags.split()]] = 1\n\n    df_tags = pd.DataFrame(tags_array)\n\n    D = {}\n    for e, i in enumerate(df_tags.columns):\n        D[i] = 'tag_' + all_tags[e]\n\n    df_tags_ = df_tags.rename(columns= D)\n    df_tags_['content_id'] = df['content_id']\n    df = df.merge(df_tags_, on='content_id')\n    \n    return df\n\ntrain_reg = create_tag_df(train_reg)\nval_reg = create_tag_df(val_reg)\ntest_reg = create_tag_df(test_reg)\n\ntrain_reg.shape, test_reg.shape, val_reg.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:02.261603Z","iopub.execute_input":"2024-04-28T12:00:02.261953Z","iopub.status.idle":"2024-04-28T12:00:04.322866Z","shell.execute_reply.started":"2024-04-28T12:00:02.261920Z","shell.execute_reply":"2024-04-28T12:00:04.321506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_reg.head()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:04.326256Z","iopub.execute_input":"2024-04-28T12:00:04.326793Z","iopub.status.idle":"2024-04-28T12:00:04.357605Z","shell.execute_reply.started":"2024-04-28T12:00:04.326749Z","shell.execute_reply":"2024-04-28T12:00:04.356393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_reg = train_reg[train_reg.bundle_id!= -10 ]\n# val_reg = val_reg[val_reg.bundle_id!= -10 ]\n# test_reg = test_reg[test_reg.bundle_id!= -10 ]\n\ntrain_reg['part'] = train_reg['part'].astype('category')\nval_reg['part'] = val_reg['part'].astype('category')\ntest_reg['part'] = test_reg['part'].astype('category')\n\ntrain_reg['bundle_id'] = train_reg['bundle_id'].astype('category')\nval_reg['bundle_id'] = val_reg['bundle_id'].astype('category')\ntest_reg['bundle_id'] = test_reg['bundle_id'].astype('category')\n\ntrain_reg.shape, test_reg.shape, val_reg.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:04.359378Z","iopub.execute_input":"2024-04-28T12:00:04.359780Z","iopub.status.idle":"2024-04-28T12:00:04.382630Z","shell.execute_reply.started":"2024-04-28T12:00:04.359746Z","shell.execute_reply":"2024-04-28T12:00:04.381162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# small sample data with train record <50 rows each item\ntrain_reg_small = train_reg[train_reg['n']<50]\ntrain_reg['n'].max()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:04.384660Z","iopub.execute_input":"2024-04-28T12:00:04.385580Z","iopub.status.idle":"2024-04-28T12:00:04.401111Z","shell.execute_reply.started":"2024-04-28T12:00:04.385526Z","shell.execute_reply":"2024-04-28T12:00:04.399661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### weighted least square\n\nfrom sklearn import linear_model, svm\nfrom sklearn.linear_model import LinearRegression\nimport matplotlib.pyplot as plt\nimport sklearn.metrics as metrics\n\ntags_col_names = train_reg.columns.tolist()[10:]\n### weighted ls regression\ntrain_reg_small = train_reg[train_reg['n']<50] #or <50, 70305 is the max n for full data\nval_reg_small = val_reg[val_reg['content_id'].isin(train_reg_small['content_id'].tolist())]\ntest_reg_small = test_reg[test_reg['content_id'].isin(train_reg_small['content_id'].tolist())]\n\nX_train = train_reg_small[['bundle_id', 'part', 'item_diff']]#+ tags_col_names]\ny_train =train_reg_small[['time_q95']]\nX_val = val_reg_small[['bundle_id', 'part', 'item_diff']]#+ tags_col_names]\ny_val =val_reg_small[['time_q95']]\nX_test = test_reg_small[['bundle_id', 'part', 'item_diff']]#+ tags_col_names]\ny_test = test_reg_small[['time_q95']]\n\nregr = LinearRegression()\nregr.fit(X_train, np.log(y_train))\n\ny_pred = regr.predict(X_train)\nresidual = (y_train - np.exp(y_pred))\nresidual = np.array(residual['time_q95'].tolist())\n\nregr = LinearRegression()\nregr.fit(X_train, np.log(y_train), 1/np.sqrt(residual**2))\ny_pred = regr.predict(X_test)\n\nexplained_variance=metrics.explained_variance_score(y_test, np.exp(y_pred))\nmean_absolute_error=metrics.mean_absolute_error(y_test, np.exp(y_pred)) \nmse=metrics.mean_squared_error(y_test, np.exp(y_pred)) \nmean_squared_log_error=metrics.mean_squared_log_error(y_test, np.exp(y_pred))\nmedian_absolute_error=metrics.median_absolute_error(y_test, np.exp(y_pred))\nr2=metrics.r2_score(y_test, y_pred)\n\nprint('explained_variance: ', round(explained_variance,4))    \nprint('mean_squared_log_error: ', round(mean_squared_log_error,4))\nprint('r2: ', round(r2,4))\nprint('MAE: ', round(mean_absolute_error,4))\nprint('MSE: ', round(mse,4))\nprint('RMSE: ', round(np.sqrt(mse),4))","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:04.403474Z","iopub.execute_input":"2024-04-28T12:00:04.405042Z","iopub.status.idle":"2024-04-28T12:00:06.423943Z","shell.execute_reply.started":"2024-04-28T12:00:04.404981Z","shell.execute_reply":"2024-04-28T12:00:06.418880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_reg_small['time_pred'] = np.exp(y_pred)\n\nratio = test_reg_small['time_pred'] / test_reg_small['time_q95'] #final_time_est\n\nra = []\nfor i in ratio.tolist():\n    if i > 2:\n        ra.append('>2')\n    elif i <= 2 and i > 1.5:\n        ra.append('1.5-2')\n    elif i <= 1.5 and i > 1.2:\n        ra.append('1.2-1.5')\n    elif i <= 1.2 and i > 1:\n        ra.append('1-1.2')\n    elif i <= 1 and i > 0.8:\n        ra.append('0.8-1')\n    elif i <= 0.8 and i > 0.5:\n        ra.append('0.5-0.8')\n    else:\n        ra.append('0-0.5')\n\ntest_reg_small['ra'] = ra\nprint(test_reg_small['ra'].value_counts())\n\n''' regression full \nra\n0.8-1      1495\n1-1.2      1357\n1.2-1.5     944\n0.5-0.8     657\n1.5-2       205\n>2           43\n0-0.5        11\nName: count, dtype: int64\n'''\n\n#regression small \n# ra\n# 0.8-1      1149\n# 1-1.2      1147\n# 1.2-1.5     840\n# 0.5-0.8     544\n# 1.5-2       203\n# >2           22\n# 0-0.5         9","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:06.434995Z","iopub.execute_input":"2024-04-28T12:00:06.436929Z","iopub.status.idle":"2024-04-28T12:00:06.480925Z","shell.execute_reply.started":"2024-04-28T12:00:06.436837Z","shell.execute_reply":"2024-04-28T12:00:06.479654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.std(np.abs(test_reg_small['time_pred'] - test_reg_small['time_q95']))","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:06.483241Z","iopub.execute_input":"2024-04-28T12:00:06.484204Z","iopub.status.idle":"2024-04-28T12:00:06.509660Z","shell.execute_reply.started":"2024-04-28T12:00:06.484154Z","shell.execute_reply":"2024-04-28T12:00:06.508431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.std(np.abs(test_reg_small['time_pred'] - test_reg_small['time_q95']))","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:06.511838Z","iopub.execute_input":"2024-04-28T12:00:06.513120Z","iopub.status.idle":"2024-04-28T12:00:06.525840Z","shell.execute_reply.started":"2024-04-28T12:00:06.513062Z","shell.execute_reply":"2024-04-28T12:00:06.524401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"X_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:06.527022Z","iopub.execute_input":"2024-04-28T12:00:06.527401Z","iopub.status.idle":"2024-04-28T12:00:06.539011Z","shell.execute_reply.started":"2024-04-28T12:00:06.527370Z","shell.execute_reply":"2024-04-28T12:00:06.537838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## xgboost\n\nfrom sklearn.model_selection import cross_val_score, PredefinedSplit\nfrom sklearn.model_selection import RepeatedKFold, GridSearchCV\nfrom sklearn.ensemble import GradientBoostingRegressor\n\n# define model\nmodel = GradientBoostingRegressor()\n\nparameters = {\n    'max_depth': [2,3,4,5,6,7,8,9,10],\n    'n_estimators': [30, 60, 100, 120],\n    'learning_rate': [0.01, 0.03, 0.05,0.07, 0.09]\n}\n\n\n#cv = RepeatedKFold(n_splits=10, n_repeats=3, random_state=1)\nps = PredefinedSplit(np.concatenate([np.zeros(X_train.shape[0])-1, np.zeros(X_val.shape[0])]))\ngrid_search = GridSearchCV(\n    estimator=model,\n    param_grid=parameters,\n    scoring = 'neg_mean_absolute_error',\n    n_jobs = 10,\n    cv = ps,\n    verbose=1\n)\n\ngrid_search.fit(pd.concat([X_train, X_val]), np.log(np.array(pd.concat([y_train['time_q95'], y_val['time_q95']]))))\n\ngrid_search.best_estimator_\n\n\n#n_scores = cross_val_score(model, X_train, y_train, scoring='neg_mean_absolute_error', cv=cv, n_jobs=-1)\n#print('MAE: %.3f (%.3f)' % (np.mean(n_scores), np.std(n_scores)))","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:06.541412Z","iopub.execute_input":"2024-04-28T12:00:06.542537Z","iopub.status.idle":"2024-04-28T12:00:43.973648Z","shell.execute_reply.started":"2024-04-28T12:00:06.542492Z","shell.execute_reply":"2024-04-28T12:00:43.971233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = GradientBoostingRegressor(max_depth=7, n_estimators=120, learning_rate=0.07)\n\nmodel.fit(X_train,np.array(y_train['time_q95']))\ny_pred_gbr = model.predict(X_test)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:43.976419Z","iopub.execute_input":"2024-04-28T12:00:43.977019Z","iopub.status.idle":"2024-04-28T12:00:44.731604Z","shell.execute_reply.started":"2024-04-28T12:00:43.976974Z","shell.execute_reply":"2024-04-28T12:00:44.729880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_reg_small['time_pred'] = y_pred_gbr\n\nratio = test_reg_small['time_pred'] / test_reg_small['time_q95'] #final_time_est\n\nra = []\nfor i in ratio.tolist():\n    if i > 2:\n        ra.append('>2')\n    elif i <= 2 and i > 1.5:\n        ra.append('1.5-2')\n    elif i <= 1.5 and i > 1.2:\n        ra.append('1.2-1.5')\n    elif i <= 1.2 and i > 1:\n        ra.append('1-1.2')\n    elif i <= 1 and i > 0.8:\n        ra.append('0.8-1')\n    elif i <= 0.8 and i > 0.5:\n        ra.append('0.5-0.8')\n    else:\n        ra.append('0-0.5')\n\ntest_reg_small['ra'] = ra\nprint(test_reg_small['ra'].value_counts())\n\n''' full\nra\n1-1.2      1725\n0.8-1      1469\n1.2-1.5     893\n0.5-0.8     401\n1.5-2       203\n>2           13\n0-0.5         8\nName: count, dtype: int64\n''' \n#small\n# ra\n# 1-1.2      1322\n# 0.8-1      1156\n# 1.2-1.5     799\n# 0.5-0.8     368\n# 1.5-2       241\n# >2           22\n# 0-0.5         6","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:44.733325Z","iopub.execute_input":"2024-04-28T12:00:44.733735Z","iopub.status.idle":"2024-04-28T12:00:44.762276Z","shell.execute_reply.started":"2024-04-28T12:00:44.733700Z","shell.execute_reply":"2024-04-28T12:00:44.760630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame({'y_red_reg': np.exp(y_pred).reshape(-1), 'y_red_gbr': y_pred_gbr}).to_csv('gbr_reg_small_pred_notag2.csv')","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:44.769792Z","iopub.execute_input":"2024-04-28T12:00:44.770283Z","iopub.status.idle":"2024-04-28T12:00:44.812324Z","shell.execute_reply.started":"2024-04-28T12:00:44.770247Z","shell.execute_reply":"2024-04-28T12:00:44.810423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.std(np.abs(test_reg_small['time_pred'] - test_reg_small['time_q95']))","metadata":{"execution":{"iopub.status.busy":"2024-04-28T12:00:44.814726Z","iopub.execute_input":"2024-04-28T12:00:44.815430Z","iopub.status.idle":"2024-04-28T12:00:44.829796Z","shell.execute_reply.started":"2024-04-28T12:00:44.815371Z","shell.execute_reply":"2024-04-28T12:00:44.828122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"in Rstudio, run the code of EBLUP\n```\nlibrary('lme4')\nlibrary('sae')\nlibrary('dplyr')\n\ndf_test_ = read.csv('test_first_half.csv')\ndf_val_ = read.csv('val_first_half.csv')\ndf_train_ = read.csv('train_first_half.csv')\ntrain_reg = read.csv('train_reg_first_half.csv')\ntest_reg = read.csv('test_reg_first_half.csv')\nval_reg = read.csv('val_reg_first_half.csv')\n\nPopn1 = read.csv('Popn1.csv')\ntime_95 = read.csv('time_95_first_half.csv')\nXmean1 = read.csv('Xmean1_.csv')\n\ndf_val_ = val_reg\ndf_test_ = test_reg\n\n###### EBLUP\ndf_train_$prior_question_elapsed_time_shift_log = log(df_train_$prior_question_elapsed_time_shift+0.01)\nelbupResults = eblupBHF(prior_question_elapsed_time_shift_log~ answered_correctly_x, \n                        dom=content_id, meanxpop=data.frame(Xmean1$content_id, Xmean1$item_diff), popnsize=Popn1, data=df_train_)\n\nresults_df = df_train_ %>% group_by(content_id) %>% summarize(Mean = mean(prior_question_elapsed_time_shift, na.rm=TRUE),\n                                                       Quant95 = quantile(prior_question_elapsed_time_shift, c(0.95)),\n                                                       sd = sd(prior_question_elapsed_time_shift, na.rm=TRUE), \n                                                       n = n(),\n                                                       mse = sd(prior_question_elapsed_time_shift, na.rm=TRUE) / sqrt(n()))\n\nresults_eblup  = data.frame (content_id  = unique(df_train_$content_id),\n                             eblup = exp(elbupResults$eblup$eblup)\n)\n\n########## merge and validate\njointdataset <- merge(results_df, results_eblup, by = 'content_id')\n#val_time_95 = df_val_ %>% group_by(content_id) %>% summarize(Quant95 = quantile(prior_question_elapsed_time_shift, c(0.95)))\n\nfreqs = NULL\nfor(z in  seq(1.5, 2, by=0.01)){\n  final_time_est = jointdataset$eblup + z*jointdataset$sd #j #1.62*jointdataset$mse*sqrt(jointdataset$n) #\n  \n  jointdataset$final_time_est = final_time_est\n  results_full = merge(jointdataset, df_val_, by = 'content_id')\n  \n  results_full <- na.omit(results_full)\n  #results_full = results_full[results_full$n>50,]\n  \n  ratio_95 = results_full['Quant95'] / results_full['time_q95'] #final_time_est\n  ratio_re = results_full['final_time_est'] / results_full['time_q95'] #final_time_est\n  \n  ra_95 = NULL\n  for(i in ratio_95$Quant95){\n    if(i>2){\n      ra_95 = c(ra_95, '>2')\n    }else if(i <= 2 && i > 1.5){\n      ra_95 = c(ra_95, '1.5-2')\n    }else if(i <= 1.5 && i > 1.2){\n      ra_95 = c(ra_95, '1.2-1.5')\n    }else if(i <= 1.2 && i > 1){\n      ra_95 = c(ra_95, '1-1.2')\n    }else if(i <= 1 && i > 0.8){\n      ra_95 = c(ra_95, '0.8-1')\n    }else if(i <= 0.8 && i > 0.5){\n      ra_95 = c(ra_95, '0.5-0.8')\n    }else{\n      ra_95 = c(ra_95, '0-0.5')\n    }\n  }\n  \n  ra_re = NULL\n  for(i in ratio_re$final_time_est){\n    \n    if(i>2){\n      ra_re = c(ra_re, '>2')\n    }else if(i <= 2 && i > 1.5){\n      ra_re = c(ra_re, '1.5-2')\n    }else if(i <= 1.5 && i > 1.2){\n      ra_re = c(ra_re, '1.2-1.5')\n    }else if(i <= 1.2 && i > 1){\n      ra_re = c(ra_re, '1-1.2')\n    }else if(i <= 1 && i > 0.8){\n      ra_re = c(ra_re, '0.8-1')\n    }else if(i <= 0.8 && i > 0.5){\n      ra_re = c(ra_re, '0.5-0.8')\n    }else{\n      ra_re = c(ra_re, '0-0.5')\n    }\n  }\n  \n  freq = as.data.frame(table(ra_re))$Freq \n  \n  freqs = cbind(freqs,  freq)\n  \n}\n\nresults = cbind(as.data.frame(table(ra_95)), freqs)\n\nre_freq_names = NULL\nfor(n in seq(1.5, 2, by=0.01)){\n  re_freq_names = c(re_freq_names, paste0('freq', n))\n}\n\ncolnames(results) = c('ra_95', 'Freq', re_freq_names)\n\n############## test\n#test_time_95 = df_test_ %>% group_by(content_id) %>% summarize(Quant95 = quantile(prior_question_elapsed_time_shift, c(0.95), na.rm=TRUE))\n##select(df_test_, content_id, time_q95) #\ntest_time_95 = df_test_[,c('content_id', 'time_q95')]\nfinal_time_est = jointdataset$eblup + 1.78*jointdataset$sd #j #1.62*jointdataset$mse*sqrt(jointdataset$n) #\n\njointdataset$final_time_est = final_time_est\nresults_full = merge(jointdataset, test_time_95, by = 'content_id')\n\n#results_full <- na.omit(results_full)\nresults_full = results_full[results_full$n<50,]\n\nratio_95 = results_full['Quant95'] / results_full['time_q95'] #final_time_est\nratio_re = results_full['final_time_est'] / results_full['time_q95'] #final_time_est\n\nra_95 = NULL\nfor(i in ratio_95$Quant95){\n  if(i>2){\n    ra_95 = c(ra_95, '>2')\n  }else if(i <= 2 && i > 1.5){\n    ra_95 = c(ra_95, '1.5-2')\n  }else if(i <= 1.5 && i > 1.2){\n    ra_95 = c(ra_95, '1.2-1.5')\n  }else if(i <= 1.2 && i > 1){\n    ra_95 = c(ra_95, '1-1.2')\n  }else if(i <= 1 && i > 0.8){\n    ra_95 = c(ra_95, '0.8-1')\n  }else if(i <= 0.8 && i > 0.5){\n    ra_95 = c(ra_95, '0.5-0.8')\n  }else{\n    ra_95 = c(ra_95, '0-0.5')\n  }\n}\n\nra_re = NULL\nfor(i in abs(ratio_re$final_time_est)){\n  \n  if(i>2){\n    ra_re = c(ra_re, '>2')\n  }else if(i <= 2 && i > 1.5){\n    ra_re = c(ra_re, '1.5-2')\n  }else if(i <= 1.5 && i > 1.2){\n    ra_re = c(ra_re, '1.2-1.5')\n  }else if(i <= 1.2 && i > 1){\n    ra_re = c(ra_re, '1-1.2')\n  }else if(i <= 1 && i > 0.8){\n    ra_re = c(ra_re, '0.8-1')\n  }else if(i <= 0.8 && i > 0.5){\n    ra_re = c(ra_re, '0.5-0.8')\n  }else{\n    ra_re = c(ra_re, '0-0.5')\n  }\n}\n\ntable(ra_re)\ntable(ra_95)\n\n```","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}