{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":50160,"databundleVersionId":7602123,"sourceType":"competition"}],"dockerImageVersionId":30646,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nimport warnings\nwarnings.simplefilter(\"ignore\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-02-25T08:57:41.852010Z","iopub.execute_input":"2024-02-25T08:57:41.852379Z","iopub.status.idle":"2024-02-25T08:57:41.857530Z","shell.execute_reply.started":"2024-02-25T08:57:41.852348Z","shell.execute_reply":"2024-02-25T08:57:41.856343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"0.85*6","metadata":{"execution":{"iopub.status.busy":"2024-02-25T08:57:41.995557Z","iopub.execute_input":"2024-02-25T08:57:41.995981Z","iopub.status.idle":"2024-02-25T08:57:42.003591Z","shell.execute_reply.started":"2024-02-25T08:57:41.995954Z","shell.execute_reply":"2024-02-25T08:57:42.002009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N_WEEKS_TOTAL = 9","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:03:22.483739Z","iopub.execute_input":"2024-02-25T10:03:22.484349Z","iopub.status.idle":"2024-02-25T10:03:22.488392Z","shell.execute_reply.started":"2024-02-25T10:03:22.484309Z","shell.execute_reply":"2024-02-25T10:03:22.487045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_normal_dist(scale, mean=0.7):\n    np.random.seed(0)\n    ary = np.random.normal(size=N_WEEKS_TOTAL, scale=scale)\n    ary = (ary - ary.mean()) / ary.std() * scale\n    ary += mean\n    return ary\n\n# test\nary = create_normal_dist(scale=0.05)\nprint(f\"[create_normal_dist] mean: {ary.mean()}, std: {ary.std()} ({ary.tolist()})\" )","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:25:07.316444Z","iopub.execute_input":"2024-02-25T10:25:07.316795Z","iopub.status.idle":"2024-02-25T10:25:07.324720Z","shell.execute_reply.started":"2024-02-25T10:25:07.316769Z","shell.execute_reply":"2024-02-25T10:25:07.323160Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef calc_high_low(scale, week, mean):\n    low = mean - scale * np.sqrt(week / (N_WEEKS_TOTAL - week))\n    high = scale * np.sqrt((N_WEEKS_TOTAL - week) / week) + mean\n    return high, low\n\ndef create_domain_shift_randomly(scale, week=6, mean=0.7):\n    high, low = calc_high_low(scale, week, mean)\n    ary = np.array([high] * week + [low] * (N_WEEKS_TOTAL - week))\n    np.random.shuffle(ary)\n    return ary\n\ndef create_domain_shift_after_nw(scale, week=6, mean=0.7):\n    high, low = calc_high_low(scale, week, mean)\n    return np.array([high] * week + [low] * (N_WEEKS_TOTAL - week))\n\n# test\nary = create_domain_shift_randomly(scale=0.05, week=6)\nprint(f\"[create_domain_shift_randomly] mean: {ary.mean()}, std: {ary.std()} ({ary.tolist()})\" )\nary = create_domain_shift_after_nw(scale=0.01, week=8)\nprint(f\"[create_domain_shift_after_nw] mean: {ary.mean()}, std: {ary.std()} ({ary.tolist()})\")","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:29:00.802754Z","iopub.execute_input":"2024-02-25T10:29:00.803115Z","iopub.status.idle":"2024-02-25T10:29:00.812331Z","shell.execute_reply.started":"2024-02-25T10:29:00.803085Z","shell.execute_reply":"2024-02-25T10:29:00.811461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_step(scale):\n    \n    residual = (np.arange(N_WEEKS_TOTAL//2+1) ** 2).sum() * 2\n    return np.sqrt(N_WEEKS_TOTAL) * scale / np.sqrt(residual)\n\ndef create_step_ary(scale, updown):\n    step = calc_step(scale=scale)\n    if updown == \"up\":\n        return np.arange(0.7-step*4, 0.7+step*4.001, step)\n    elif updown == \"down\":\n        return np.arange(0.7+step*4, 0.7-step*4.001, -step)\n    \n# test\nary = create_step_ary(0.05, updown=\"up\")\nprint(f\"[create_step_ary] mean: {ary.mean()}, std: {ary.std()} ({ary.tolist()})\")","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:29:00.954719Z","iopub.execute_input":"2024-02-25T10:29:00.955086Z","iopub.status.idle":"2024-02-25T10:29:00.964241Z","shell.execute_reply.started":"2024-02-25T10:29:00.955059Z","shell.execute_reply":"2024-02-25T10:29:00.962940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.random.seed(0)\n\n# total: 6.3 (0.7 * 9weeks)\ngini_in_time_dict = {\n    \"constant value\": np.array([0.7] * 9),\n}\n\nfor std in [0.01, 0.025, 0.05]:\n    gini_in_time_dict[f\"monotonically decreasing(std={std})\"] = create_step_ary(std, updown=\"down\")\n    gini_in_time_dict[f\"monotonically increasing(std={std})\"] = create_step_ary(std, updown=\"up\")\n    gini_in_time_dict[f\"normal_dist(std={std})\"] = create_normal_dist(scale=std)\n    for week in [6, 8]:\n        gini_in_time_dict[f\"domain_shift_after_{week}w(std={std})\"] = create_domain_shift_after_nw(scale=std, week=week)\n        gini_in_time_dict[f\"domain_shift_randomly_{week}w(std={std})\"] = create_domain_shift_randomly(scale=std, week=week)\n    \ngini_in_time_dict","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:32:41.945801Z","iopub.execute_input":"2024-02-25T10:32:41.946146Z","iopub.status.idle":"2024-02-25T10:32:41.961380Z","shell.execute_reply.started":"2024-02-25T10:32:41.946117Z","shell.execute_reply":"2024-02-25T10:32:41.960571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 4))\nfor k, v in gini_in_time_dict.items():\n    assert round(v.sum(), 1) == 6.3, f\"{k}: {v.sum()}, {v}\"\n    plt.plot(v, label=k)\nplt.ylim(0.5, 0.9)\nplt.xlabel(\"time\")\nplt.ylabel(\"score\")\nplt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:32:42.103727Z","iopub.execute_input":"2024-02-25T10:32:42.104087Z","iopub.status.idle":"2024-02-25T10:32:42.477697Z","shell.execute_reply.started":"2024-02-25T10:32:42.104057Z","shell.execute_reply":"2024-02-25T10:32:42.476319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def gini_stability(gini_in_time, w_fallingrate=88.0, w_resstd=-0.5):\n    x = np.arange(len(gini_in_time))\n    y = gini_in_time\n    a, b = np.polyfit(x, y, 1)\n    y_hat = a*x + b\n    residuals = y - y_hat\n    res_std = np.std(residuals)\n    avg_gini = np.mean(gini_in_time)\n    return avg_gini + w_fallingrate * min(0, a) + w_resstd * res_std\n\ndef metric_jacobyjaeger(gini_in_time, exponent=69, ma_len=4):\n    x = np.cumsum(gini_in_time, axis=0)\n    x = np.concatenate([0*x[:1], x], 0)\n    scores = -np.mean(-np.log(np.maximum((x[ma_len:] - x[:-ma_len])/ma_len, 1e-5))**exponent)\n    return scores \n\ndef metric_seifachour12(gini_in_time, alpha=1028, beta=1.97, gamma=0.49):\n    return (np.abs(beta*(1/beta*np.mean(gini_in_time) + (gamma-np.std(gini_in_time)))-1))**(1/alpha) \n\ndef metric_kononenko(gini_in_time, c=0.90):\n    return np.mean(gini_in_time) - c*np.std(gini_in_time)\n\ndef metric_davutpolat(gini_in_time, c=1.82):\n    max_gini = gini_in_time[0]\n    cost = 0\n    for i in range(1, len(gini_in_time)):\n        max_gini = max(max_gini, gini_in_time[i])\n        cost += max(0, max_gini - gini_in_time[i])**2\n    # return np.mean(gini_in_time) + c/(len(gini_in_time)-1)*cost   \n    return np.mean(gini_in_time) - c/(len(gini_in_time)-1)*cost   \n\ndef gini_at7459(gini_in_time, w_fallingrate=88.0, w_resstd=-0.5, f=8):\n    w_fallingrate /= f + 1\n\n    x = np.arange(len(gini_in_time))\n    y = gini_in_time\n    a, b = np.polyfit(x, y, 1)\n    y_hat = a*x + b\n    residuals = y - y_hat\n    res_std = np.std(residuals)\n    avg_gini = np.mean(gini_in_time)\n    far = a\n    for s in range(1,f):\n        start_index = len(x) // f * (s)\n        end_index = len(x) // f * (s+1)\n        x_second_fifth = x[start_index:end_index]\n        y_second_fifth = gini_in_time[start_index:end_index]\n        x1 = x[start_index:end_index]\n        y1 = gini_in_time[start_index:end_index]\n        a1, b1 = np.polyfit(x1, y1, 1)\n        far += min(0,a1)\n\n    return avg_gini + w_fallingrate * (far) + w_resstd * res_std \n\ndef metric_mean(gini_in_time):\n    return np.mean(gini_in_time)  \n\ndef eval_metric(gini_in_time):\n    ret = {}\n    for func in [\n        gini_stability,\n        metric_jacobyjaeger,\n        metric_seifachour12,\n        metric_kononenko,\n        metric_davutpolat,\n        gini_at7459,\n        metric_mean,\n    ]:\n        ret[func.__name__] = func(gini_in_time)\n    return ret\n        \n","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:32:42.479823Z","iopub.execute_input":"2024-02-25T10:32:42.480199Z","iopub.status.idle":"2024-02-25T10:32:42.500615Z","shell.execute_reply.started":"2024-02-25T10:32:42.480169Z","shell.execute_reply":"2024-02-25T10:32:42.499248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_result = {}\nfor name, gini_in_time in gini_in_time_dict.items():\n    df_result[name] = eval_metric(gini_in_time)","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:32:42.542549Z","iopub.execute_input":"2024-02-25T10:32:42.542918Z","iopub.status.idle":"2024-02-25T10:32:42.568553Z","shell.execute_reply.started":"2024-02-25T10:32:42.542891Z","shell.execute_reply":"2024-02-25T10:32:42.567562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_result = pd.DataFrame(df_result).T","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:32:42.684032Z","iopub.execute_input":"2024-02-25T10:32:42.684977Z","iopub.status.idle":"2024-02-25T10:32:42.690330Z","shell.execute_reply.started":"2024-02-25T10:32:42.684942Z","shell.execute_reply":"2024-02-25T10:32:42.689302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_result","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:32:43.034979Z","iopub.execute_input":"2024-02-25T10:32:43.035896Z","iopub.status.idle":"2024-02-25T10:32:43.055704Z","shell.execute_reply.started":"2024-02-25T10:32:43.035862Z","shell.execute_reply":"2024-02-25T10:32:43.054602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Consider about prefer model\nIf the sum and standard deviation of scores and are the same, the model would consider the following order preferable for business:\n\n1. Normal distribution\n2. Monotonically increasing or monotonically decreasing\n3. Domain shift","metadata":{}},{"cell_type":"code","source":"for i, row in df_result.T.iterrows():\n    print(\"--------------------\")\n    print(row.name)\n    print(\"--------------------\")\n    row = row.sort_values(ascending=False)\n    display(row)","metadata":{"execution":{"iopub.status.busy":"2024-02-25T10:42:42.670668Z","iopub.execute_input":"2024-02-25T10:42:42.671086Z","iopub.status.idle":"2024-02-25T10:42:42.710705Z","shell.execute_reply.started":"2024-02-25T10:42:42.671051Z","shell.execute_reply":"2024-02-25T10:42:42.710058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}