{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.11.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":101849,"databundleVersionId":12846694,"sourceType":"competition"},{"sourceId":12444802,"sourceType":"datasetVersion","datasetId":7850239}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Competition metric in detail\n\n## Summary\n* We can use the official metric code with `fgs_weight = 0.4 / 1.95 * 282`.\n* Note that L_ref is *not* a constant, it depends on y_true, calculated for each n_planets × wavelenth, and has shape (n_planets, 283), according to the official code.\n\n\n## Competition metric\n\nEvaluation section: https://www.kaggle.com/competitions/ariel-data-challenge-2025/overview/evaluation\n\nGaussian log likelihood:\n$$\\mathcal{L}(y_\\mathrm{true}, \\mu_\\mathrm{pred}, \\sigma_\\mathrm{pred})  = -\\frac{1}{2} \\left[　\\log(2\\pi \\sigma_\\mathrm{pred}^2) + \\frac{ \\left( y_\\mathrm{true} - \\mu_\\mathrm{pred} \\right)^2 }{ \\sigma_\\mathrm{pred}^2 } \\right]$$\n\nScores for n_planets × wavelength:\n\n$$\\mathrm{scores}(y_\\mathrm{true}, \\mu_\\mathrm{pred}, \\sigma_\\mathrm{pred}) = \\frac{\\mathcal{L}(y_\\mathrm{true}, \\mu_\\mathrm{pred}, \\sigma_\\mathrm{pred}) - \\mathcal{L}_\\mathrm{ref}(y_\\mathrm{true})} {\\mathcal{L}_\\mathrm{ideal} - \\mathcal{L}_\\mathrm{ref}(y_\\mathrm{true})} $$\n\nThe score is the weighted average of the scores\n\n> FGS1: 2 × (0.2/1) = 0.4<br/>\n> AIRS-Ch0: 1.95/282 ≈ 0.0069 per spectral point\n\n\nOfficial metric code:\n\nhttps://www.kaggle.com/code/metric/ariel-gaussian-log-likelihood\n\nSee also a great discussion by AmbrosM for Ariel 2024:\n\nhttps://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/528114\n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport sys\nimport scipy\n\nsys.path.append('/kaggle/input/ariel2025-public/score')\nimport official_metric\n\nINPUT_DIR = '/kaggle/input/ariel-data-challenge-2025/'","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.400177Z","iopub.execute_input":"2025-07-12T11:50:05.400404Z","iopub.status.idle":"2025-07-12T11:50:05.416829Z","shell.execute_reply.started":"2025-07-12T11:50:05.400385Z","shell.execute_reply":"2025-07-12T11:50:05.415958Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Use official metric","metadata":{}},{"cell_type":"code","source":"# Label y_true for 1100 training data\ntrain = pd.read_csv(INPUT_DIR + '/train.csv')\ny_true = train.iloc[:, 1:].to_numpy()\n\n# Submission file for 1100 training data\nsubmission = pd.read_csv('/kaggle/input/ariel2025-public/score/submission_train.csv')\nmu_pred = submission.iloc[:, 1:284].to_numpy()\nsigma_pred = submission.iloc[:, 284:].to_numpy()\n\n# Assume planet_id in same order\nassert (train['planet_id'] == submission['planet_id']).all()\n\ny_true.shape, mu_pred.shape, sigma_pred.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.417756Z","iopub.execute_input":"2025-07-12T11:50:05.417993Z","iopub.status.idle":"2025-07-12T11:50:05.762036Z","shell.execute_reply.started":"2025-07-12T11:50:05.417973Z","shell.execute_reply":"2025-07-12T11:50:05.761317Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"naive_mean = np.mean(y_true)\nnaive_sigma = np.std(y_true, ddof=1)  # ddof=1 makes negligible difference but I argue that this is statistally correct\n\nprint('naive_mean:  %.8e' % naive_mean)\nprint('naive_sigma: %.8e' % naive_sigma)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.762917Z","iopub.execute_input":"2025-07-12T11:50:05.763244Z","iopub.status.idle":"2025-07-12T11:50:05.772523Z","shell.execute_reply.started":"2025-07-12T11:50:05.763216Z","shell.execute_reply":"2025-07-12T11:50:05.771465Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# airs_weight=1 is fixed in the official metric, set fgs_weight relative to this airs_weight\nfgs_weight = 0.4 / 1.95 * 282\nprint('fgs_weight: %.8e' % fgs_weight)\n\ns = official_metric.score(train.copy(), submission.copy(), 'planet_id', naive_mean, naive_sigma, fgs_weight=fgs_weight)\n# Because score function drop planet_id column, I pass the copies; in case you use them later.\n\nprint('Local CV:  %.4f' % s)\nprint('Public LB: %.3f' % 0.303)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.773622Z","iopub.execute_input":"2025-07-12T11:50:05.773876Z","iopub.status.idle":"2025-07-12T11:50:05.903246Z","shell.execute_reply.started":"2025-07-12T11:50:05.773853Z","shell.execute_reply":"2025-07-12T11:50:05.902217Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Because training set is not exactly same as the hidden test set, Local CV and Public LB are not exactly the same.","metadata":{}},{"cell_type":"markdown","source":"# Check score\n\nI made several submissions with same `mu_pred` but various `sigma_pred` to confirm the competition metric and the value of L_ref.\n\n**You don't have to read below. You can just use the official metric code.** The following is for those who doubt and want to confirm exactly what the competition metric is. One thing I misunderstood (and I guess the code is doing what the organizer does not expect)\nis that L_ref is not a constant and depends on planet_id × wavelength.\n\n1. I thought L_ref was a constant and reverse engineered the L_ref value from Public LB.\n2. This effective constant L_ref = 3.23542514 explains my Public LB scores.\n3. But it is not equal to the mean L_ref = 3.1222 over training set.\n4. I understand the cause of the discrepancy (2) and (3):\n   - L_ref is not a constant,\n   - and, score = (L - L_ref) / (L_ideal - L_ref) is nonlinear in L_ref.\n\n* Other parameters, L_ideal and weights for FGS1 and AIRS agree with the description.\n","metadata":{}},{"cell_type":"code","source":"# My Public LB data\nsigma_fgs  = np.array([ 1e-3, 0.5e-3,  4e-3,  2e-3,  2e-3,  4e-3])\nsigma_airs = np.array([ 1e-3, 0.5e-3,  4e-3,  2e-3,  1e-3,  1e-3])\nscores     = np.array([0.303,  0.226, 0.175, 0.255, 0.300, 0.291])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.904296Z","iopub.execute_input":"2025-07-12T11:50:05.90461Z","iopub.status.idle":"2025-07-12T11:50:05.910962Z","shell.execute_reply.started":"2025-07-12T11:50:05.904577Z","shell.execute_reply":"2025-07-12T11:50:05.90991Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sigma_ideal_fgs = 1e-6\nsigma_ideal_airs = 1e-5\nsigma_ideal = np.array([sigma_ideal_fgs, ] + [sigma_ideal_airs, ] * 282).reshape(1, 283)\n\nL_ideal = -0.5 * np.log(2 * np.pi * sigma_ideal ** 2)\n\nprint('L_ideal: %.8f %.8f' % (L_ideal[0, 0], L_ideal[0, 1]), sigma_ideal.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.911985Z","iopub.execute_input":"2025-07-12T11:50:05.912307Z","iopub.status.idle":"2025-07-12T11:50:05.932003Z","shell.execute_reply.started":"2025-07-12T11:50:05.912277Z","shell.execute_reply":"2025-07-12T11:50:05.931123Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"GLL_mean = scipy.stats.norm.logpdf(y_true,\n                                   loc=naive_mean * np.ones_like(y_true),\n                                   scale=naive_sigma * np.ones_like(y_true))\nGLL_mean.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.935559Z","iopub.execute_input":"2025-07-12T11:50:05.935802Z","iopub.status.idle":"2025-07-12T11:50:05.985725Z","shell.execute_reply.started":"2025-07-12T11:50:05.935781Z","shell.execute_reply":"2025-07-12T11:50:05.984864Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"I first thought that L_ref is a constant, mean over all training data, but that is not correct.","metadata":{}},{"cell_type":"code","source":"print('L_ref not the constant: %.4f' % GLL_mean.mean())\nprint('or wavelentgh-dependent mean: %r...' % [float(x) for x in GLL_mean.mean(axis=0)[:4]])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.986499Z","iopub.execute_input":"2025-07-12T11:50:05.986723Z","iopub.status.idle":"2025-07-12T11:50:05.992948Z","shell.execute_reply.started":"2025-07-12T11:50:05.986704Z","shell.execute_reply":"2025-07-12T11:50:05.99212Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Fit public LB scores\n\nNote that I am not assuming training set is equal to test set.\nI only assume that the functional form of the metric is known and fit for unknown Mean Squared Errors (MSE) and one more parameter L_ref.\n","metadata":{}},{"cell_type":"code","source":"# Data\nsigma_fgs  = np.array([ 1e-3, 0.5e-3,  4e-3,  2e-3,  2e-3,  4e-3])\nsigma_airs = np.array([ 1e-3, 0.5e-3,  4e-3,  2e-3,  1e-3,  1e-3])\nscores     = np.array([0.303,  0.226, 0.175, 0.255, 0.300, 0.291])\n\n# Constant\nsigma_perf_fgs = 1e-6\nsigma_perf_airs = 1e-5\n\nL_perf_fgs = -0.5 * np.log(2 * np.pi * sigma_perf_fgs ** 2)\nL_perf_airs = -0.5 * np.log(2 * np.pi * sigma_perf_airs ** 2)\nprint('L_perf:', L_perf_fgs, L_perf_airs)\n\n\ndef compute_score(sigma_fgs, sigma_airs, mse_fgs, mse_airs, L_ref):\n    \"\"\"\n    Args:\n      sigma_fgs (array):  submission data (n, )\n      sigma_airs (array): submission data (n, )\n      mse_fgs, mse_airs (float): mse = [(mu_pred - mu_true) / 1e-3] ** 2\n      w_fgs (float): unnormalized weight for FGS1; w_airs = 1\n\n    Returns:\n      scores (array): (n, )\n    \"\"\"\n    w_airs = 1.95\n    w_fgs = 0.4\n    L_fgs  = -0.5 * (np.log(2 * np.pi * sigma_fgs ** 2)  + mse_fgs * (1e-3 / sigma_fgs) ** 2)\n    L_airs = -0.5 * (np.log(2 * np.pi * sigma_airs ** 2) + mse_airs * (1e-3 / sigma_airs) ** 2)\n\n    scores = (w_fgs / (L_perf_fgs - L_ref) * (L_fgs - L_ref) + w_airs / (L_perf_airs - L_ref) * (L_airs - L_ref)) / (w_fgs + w_airs)\n    return scores\n\n\ndef f_opt(x):\n    # Optimize unknown parameters\n    mse_fgs, mse_airs, L_ref = x\n    scores_theory = compute_score(sigma_fgs, sigma_airs, mse_fgs, mse_airs, L_ref)\n    loss = np.mean((scores_theory - scores) ** 2)\n    return loss\n\nopt = scipy.optimize.minimize(f_opt, x0=(1, 1, 3))\nprint('opt parameters', opt.x)\nprint('rmse', np.sqrt(opt.fun))\n\n\nmse_fgs, mse_airs, L_ref = opt.x\n\nsigma_smooth = np.linspace(5e-4, 4e-3, 101)\nones = 1e-3 * np.ones_like(sigma_smooth)\n\nscores_smooth  = compute_score(sigma_smooth, sigma_smooth, mse_fgs, mse_airs, L_ref)\nscores_smooth2 = compute_score(sigma_smooth,         ones, mse_fgs, mse_airs, L_ref)\n\nplt.title('Scores for fixed $\\\\mu_\\\\mathrm{pred}$; $\\\\mathcal{L}_\\\\mathrm{ref} = %.4f$' % L_ref)\nplt.plot(sigma_smooth, scores_smooth, alpha=0.8, label='$\\\\sigma_\\\\mathrm{FGS} = \\\\sigma_\\\\mathrm{AIRS} = \\\\sigma$')\nplt.plot(sigma_smooth, scores_smooth2, alpha=0.8, label='$\\\\sigma_\\\\mathrm{FGS} = \\\\sigma; \\\\, \\\\sigma_\\\\mathrm{AIRS}=10^{-3}$')\nplt.plot(sigma_fgs, scores, 'x', color='black', label='Local CV scores')\nplt.xlabel('$\\\\sigma$')\nplt.ylabel('Score')\nplt.legend()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:05.993903Z","iopub.execute_input":"2025-07-12T11:50:05.994223Z","iopub.status.idle":"2025-07-12T11:50:06.329703Z","shell.execute_reply.started":"2025-07-12T11:50:05.994195Z","shell.execute_reply":"2025-07-12T11:50:06.328729Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('Effective constant L_ref:', L_ref)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:06.331043Z","iopub.execute_input":"2025-07-12T11:50:06.331399Z","iopub.status.idle":"2025-07-12T11:50:06.336706Z","shell.execute_reply.started":"2025-07-12T11:50:06.331366Z","shell.execute_reply":"2025-07-12T11:50:06.335678Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"If L_ref were constant, it would be 3.233, not 3.1222\n\nI would say, the other parameters\n\n* L_ideal_fgs, L_ideal_airs, w_fgs/w_airs\n\nare correct because they can explain the public LB sigma dependence.","metadata":{}},{"cell_type":"markdown","source":"## Understand effective L_ref\nThe reason that mean $L_\\mathrm{ref}$ is not correct is clear -- it is **not** constant.\n\nScores are nonlinear function of $\\mathcal{L}_\\mathrm{ref}$,\n\n$$\n\\mathrm{scores}\n  = \\frac{\\mathcal{L} - \\mathcal{L}_\\mathrm{ref}}\n         {\\mathcal{L}_\\mathrm{ideal} - \\mathcal{L}_\\mathrm{ref}},\n$$\n\ntherefore, mean $\\mathcal{L}_\\mathrm{ref}$ does not give mean score.\n\nI am going to show that mean coefficients explain the effective L_ref = 3.23398594:\n\n$$\n\\mathrm{scores} = a \\mathcal{L} - b\n$$\n\n$$\n\\mathrm{score} \\approx \\langle a \\rangle \\langle \\mathcal{L} \\rangle - \\langle b \\rangle\n$$\n\nNote that this approximation is assuming log-likelihood L is uncorrelated with L_ref, which is not always true and depends on your prediction.","metadata":{}},{"cell_type":"code","source":"GLL_pred = scipy.stats.norm.logpdf(y_true, loc=mu_pred, scale=sigma_pred)\nGLL_ideal = scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_ideal * np.ones_like(y_true))\nGLL_ref = scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true))\n\nGLL_pred.shape, GLL_ideal.shape, GLL_ref.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:06.337771Z","iopub.execute_input":"2025-07-12T11:50:06.338116Z","iopub.status.idle":"2025-07-12T11:50:06.40423Z","shell.execute_reply.started":"2025-07-12T11:50:06.338085Z","shell.execute_reply":"2025-07-12T11:50:06.403378Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# score = a L - b\na = 1 / (GLL_ideal - GLL_ref)\nb = GLL_ref / (GLL_ideal - GLL_ref)\n\na.shape, b.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:06.405296Z","iopub.execute_input":"2025-07-12T11:50:06.405576Z","iopub.status.idle":"2025-07-12T11:50:06.413325Z","shell.execute_reply.started":"2025-07-12T11:50:06.405554Z","shell.execute_reply":"2025-07-12T11:50:06.412206Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Plot L_ref histogram\nbins = np.linspace(-15, 10, 101)\n\nplt.figure(figsize=(6, 2))\nplt.title('L_ref is not a constant')\nplt.xlabel('L_ref')\nplt.hist(GLL_ref[:, 0], bins, histtype='step', density=True, label='FGS1')\nplt.hist(GLL_ref[:, 1:].flatten(), bins, histtype='step', density=True, label='AIRS')\nplt.legend(loc=2, frameon=False)\nplt.show()\n\n# Plot a, b, where\n# score = (L - L_ref) / (L_ideal - L_ref) = a L - b\nplt.figure(figsize=(8, 2.5))\nplt.suptitle('$\\\\mathrm{score} = a \\\\mathcal{L} - b$')\nplt.subplot(1, 2, 1)\nplt.xlabel('a')\nbins = np.linspace(0.03, 0.15, 101)\nplt.hist(a[:, 0], bins, label='FGS1', density=True, histtype='step')\nplt.hist(a[:, 1:].flatten(), bins, label='AIRS', density=True, histtype='step')\nplt.legend(loc=2, frameon=False)\n\nplt.subplot(1, 2, 2)\nplt.xlabel('b')\nbins = np.linspace(-0.6, 0.6, 101)\nplt.hist(b[:, 0], bins, label='FGS1', density=True, histtype='step')\nplt.hist(b[:, 1:].flatten(), bins, label='AIRS', density=True, histtype='step')\n\nplt.tight_layout()\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:51:12.505674Z","iopub.execute_input":"2025-07-12T11:51:12.505964Z","iopub.status.idle":"2025-07-12T11:51:13.075225Z","shell.execute_reply.started":"2025-07-12T11:51:12.505945Z","shell.execute_reply":"2025-07-12T11:51:13.074148Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"planet_id × wavelength dependent L_ref means Gaussian log likelihoods are weighted differently depending on y_true; planet with small `a(y_true)` have small weight in the final score and vise versa.","metadata":{}},{"cell_type":"code","source":"print('FGS1')\nprint('  Mean      a:', np.mean(a[:, 0]))\nprint('  Effective a:', 1 / (L_ideal[0, 0] - L_ref))\nprint('')\nprint('  Mean      b:', np.mean(b[:, 0]))\nprint('  Effective b:', L_ref / (L_ideal[0, 0] - L_ref))\n\nprint('')\nprint('AIRS')\nprint('  Mean      a:', np.mean(a[:, 1:]))\nprint('  Effective a:', 1 / (L_ideal[0, 1] - L_ref))\nprint('')\nprint('  Mean      b:', np.mean(b[:, 1:]))\nprint('  Effective b:', L_ref / (L_ideal[0, 1] - L_ref))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-12T11:50:06.935723Z","iopub.execute_input":"2025-07-12T11:50:06.935968Z","iopub.status.idle":"2025-07-12T11:50:06.946288Z","shell.execute_reply.started":"2025-07-12T11:50:06.935948Z","shell.execute_reply":"2025-07-12T11:50:06.945593Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"These numbers show that constant effective L_ref=3.233 can fit the Public LB scores because it gives effective coefficients `a_eff(L_ref), b_eff(L_ref)` that are close to average coefficients `mean(a), mean(b)`.","metadata":{}},{"cell_type":"markdown","source":"# Conclusion\n\n* I think I understand the competition metric and all the paramters.\n* When you use the competition metric, use `fgs_weight = 0.4 / 1.95 * 282` in the official discription, which is not the default paramter fgs_weight=1 in the code.\n* You can use default paramters for `fsg_sigma_true=1e-6`, `airs_sigma_true=1e-5`.\n* The reason I said I understand is that I understand the change in Public LB as I change sigma_pred.\n* I understant why fitting paramter L_ref is different from mean L_ref; it is due to the non-linear mapping L_ref -> score.\n","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}