{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":81933,"databundleVersionId":9643020,"sourceType":"competition"}],"dockerImageVersionId":30786,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import polars as pl\nfrom pathlib import Path\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import Dataset, DataLoader\nimport numpy as np\nimport pandas as pd\nfrom collections import defaultdict\nfrom catboost import CatBoostRegressor, Pool\nimport lightgbm as lgb\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport numpy as np\nimport pandas as pd\nimport polars as pl\nimport sklearn.neighbors, sklearn.metrics, sklearn.preprocessing\nfrom sklearn.model_selection import train_test_split, KFold\nimport xgboost as xgb\nimport polars.selectors as cs\nfrom sklearn.metrics import cohen_kappa_score, ConfusionMatrixDisplay\nimport numpy as np\nfrom colorama import Fore, Style\nimport random\nfrom scipy.stats import mode\n\nimport polars as pl\nimport numpy as np\nfrom sklearn.preprocessing import StandardScaler\nfrom tqdm import tqdm\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport gc\nimport itertools\nimport pickle\nimport re\nimport time\nimport os\nimport logging \nfrom scipy.stats.mstats import winsorize\nfrom sklearn.model_selection import StratifiedKFold\n\n# --------- config\n\nN_THREADS = 3\nN_BAGS = 10\nN_FOLDS = 5\nN_SEEDS = 1\nSEED = 0\nGPU = 0\n\nos.environ['PYTHONHASHSEED'] = str(SEED)\nrandom.seed(SEED)\nnp.random.seed(SEED)\nos.environ['POLARS_MAX_THREADS'] = str(N_THREADS)\n\nFIX_SII = False","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:26.156669Z","iopub.execute_input":"2024-11-06T01:51:26.157775Z","iopub.status.idle":"2024-11-06T01:51:26.177523Z","shell.execute_reply.started":"2024-11-06T01:51:26.157728Z","shell.execute_reply":"2024-11-06T01:51:26.176000Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"Conformal classifiers, regressors, and predictive systems (crepes) extras\n\nFunctions for generating non-conformity scores and Mondrian categories\n(bins), and classes for generating difficulty estimates and Mondrian \ncategorizers, with and without out-of-bag predictions.\n\nAuthor: Henrik Boström (bostromh@kth.se)\n\nCopyright 2024 Henrik Boström\n\nLicense: BSD 3 clause\n\n\"\"\"\n\nimport numpy as np\nimport pandas as pd\n\nfrom sklearn.neighbors import NearestNeighbors\nfrom sklearn.preprocessing import MinMaxScaler\n\ndef hinge(X_prob, classes=None, y=None):\n \"\"\"\n Computes non-conformity scores for conformal classifiers.\n\n Parameters\n ----------\n X_prob : array-like of shape (n_samples, n_classes)\n predicted class probabilities\n classes : array-like of shape (n_classes,), default=None\n class names\n y : array-like of shape (n_samples,), default=None\n correct target values\n\n Returns\n -------\n scores : ndarray of shape (n_samples,) or (n_samples, n_classes)\n non-conformity scores. The shape is (n_samples, n_classes)\n if classes and y are None.\n\n Examples\n --------\n Assuming that ``X_prob`` is an array with predicted probabilities and\n ``classes`` and ``y`` are vectors with the class names (in order) and\n correct class labels, respectively, the non-conformity scores are generated \n by:\n\n .. code-block:: python\n\n from crepes.extras import hinge\n \n alphas = hinge(X_prob, classes, y)\n\n The above results in that ``alphas`` is assigned a vector of the same length\n as ``X_prob`` with a non-conformity score for each object, here \n defined as 1 minus the predicted probability for the correct class label.\n These scores can be used when fitting a :class:`.ConformalClassifier` or\n calibrating a :class:`.WrapClassifier`. Non-conformity scores for test \n objects, for which ``y`` is not known, can be obtained from the corresponding\n predicted probabilities (``X_prob_test``) by:\n\n .. code-block:: python\n\n alphas_test = hinge(X_prob_test)\n\n The above results in that ``alphas_test`` is assigned an array of the same\n shape as ``X_prob_test`` with non-conformity scores for each class in the \n columns for each test object.\n \"\"\"\n if y is not None:\n if isinstance(y, pd.Series):\n y = y.values\n class_indexes = np.array(\n [np.argwhere(classes == y[i])[0][0] for i in range(len(y))])\n result = 1-X_prob[np.arange(len(y)),class_indexes]\n else:\n result = 1-X_prob\n return result\n\ndef margin(X_prob, classes=None, y=None):\n \"\"\"Computes non-conformity scores for conformal classifiers.\n\n Parameters\n ----------\n X_prob : array-like of shape (n_samples, n_classes)\n predicted class probabilities\n classes : array-like of shape (n_classes,), default=None\n class names\n y : array-like of shape (n_samples,), default=None\n correct target values\n\n Returns\n -------\n scores : ndarray of shape (n_samples,) or (n_samples, n_classes)\n non-conformity scores. The shape is (n_samples, n_classes)\n if classes and y are None.\n\n Examples\n --------\n Assuming that ``X_prob`` is an array with predicted probabilities and\n ``classes`` and ``y`` are vectors with the class names (in order) and\n correct class labels, respectively, the non-conformity scores are generated \n by:\n\n .. code-block:: python\n\n from crepes.extras import margin\n \n alphas = margin(X_prob, classes, y)\n\n The above results in that ``alphas`` is assigned a vector of the same length \n as ``X_prob`` with a non-conformity score for each object, here\n defined as the highest predicted probability for a non-correct class label \n minus the predicted probability for the correct class label. These scores can\n be used when fitting a :class:`.ConformalClassifier` or calibrating a \n :class:`.WrapClassifier`. Non-conformity scores for test objects, for which \n ``y`` is not known, can be obtained from the corresponding predicted \n probabilities (``X_prob_test``) by:\n\n .. code-block:: python\n\n alphas_test = margin(X_prob_test)\n\n The above results in that ``alphas_test`` is assigned an array of the same\n shape as ``X_prob_test`` with non-conformity scores for each class in the \n columns for each test object.\n\n \"\"\"\n if y is not None:\n if isinstance(y, pd.Series):\n y = y.values\n class_indexes = np.array(\n [np.argwhere(classes == y[i])[0][0] for i in range(len(y))])\n result = np.array([\n (np.max(X_prob[i, [j != class_indexes[i]\n for j in range(X_prob.shape[1])]])\n - X_prob[i, class_indexes[i]]) for i in range(len(X_prob))])\n else:\n result = np.array([\n [(np.max(X_prob[i, [j != c for j in range(X_prob.shape[1])]])\n - X_prob[i, c]) for c in range(X_prob.shape[1])]\n for i in range(len(X_prob))])\n return result\n\ndef binning(values, bins=10):\n \"\"\"\n Provides bins for a set of values.\n\n Parameters\n ----------\n values : array-like of shape (n_samples,)\n set of values\n bins : int or array-like of shape (n_bins,), default=10\n number of bins to use for equal-sized binning or threshold values \n to use for binning\n \n Returns\n -------\n assigned_bins : array-like of shape (n_samples,)\n bins to which values have been assigned\n boundaries : array-like of shape (bins+1,)\n threshold values for the bins; the first is always -np.inf and\n the last is np.inf. Returned only if bins is an int.\n\n Examples\n --------\n Assuming that ``sigmas`` is a vector with difficulty estimates,\n then Mondrian categories (bins) can be formed by finding thresholds\n for 20 equal-sized bins by:\n\n .. code-block:: python\n\n from crepes.extras import binning\n \n bins, bin_thresholds = binning(sigmas, bins=20)\n\n The above results in that ``bins`` is assigned a vector\n of the same length as ``sigmas`` with label names (integers\n from 0 to 19), while ``bin_thresholds`` define the boundaries\n for the bins. The latter can be used to assign bin labels\n to another vector, e.g., ``sigmas_test``, by providing the thresholds \n as input to :meth:`binning`:\n\n .. code-block:: python\n \n test_bins = binning(sigmas_test, bins=bin_thresholds)\n\n Here the output is just a vector ``test_bins`` with label names\n of the same length as ``sigmas_test``.\n\n Note\n ----\n A very small random number is added to each value when forming bins\n for the purpose of tie-breaking.\n \"\"\"\n mod_values = values+np.random.rand(len(values))*1e-9\n # Adding a very small random number, which a.s. avoids ties\n # without affecting performance\n if isinstance(bins, int):\n assigned_bins, bin_boundaries = pd.qcut(mod_values,bins,\n labels=False,retbins=True,\n duplicates=\"drop\",\n precision=12)\n bin_boundaries[0] = -np.inf\n bin_boundaries[-1] = np.inf\n return assigned_bins, bin_boundaries\n else:\n assigned_bins = pd.cut(mod_values,bins,labels=False,retbins=False)\n return assigned_bins\n\nclass MondrianCategorizer():\n \"\"\"\n A MondrianCategorizer outputs categories for objects to be used by \n Mondrian conformal classifiers, regressors and predictive systems.\n \"\"\"\n \n def __init__(self):\n self.fitted = False\n self.f = None\n self.de = None\n self.learner = None\n self.oob = False\n self.bin_thresholds = None\n\n def __repr__(self):\n if self.f is not None:\n return (f\"MondrianCategorizer(fitted={self.fitted}, \"\n f\"f={self.f.__name__}, no_bins={len(self.bin_thresholds)-1})\")\n elif self.de is not None:\n return (f\"MondrianCategorizer(fitted={self.fitted}, \"\n f\"de={self.de}, no_bins={len(self.bin_thresholds)-1})\")\n elif self.learner is not None:\n return (f\"MondrianCategorizer(fitted={self.fitted}, \"\n f\"learner={self.learner}, no_bins={len(self.bin_thresholds)-1})\")\n else:\n return f\"MondrianCategorizer(fitted={self.fitted})\"\n \n def fit(self, X=None, f=None, de=None, learner=None, no_bins=10, oob=False):\n \"\"\"\n Fit Mondrian categorizer.\n\n Parameters\n ----------\n X : array-like of shape (n_samples, n_features), default=None\n set of objects\n f : function which given an array-like of shape (n_samples, n_features)\n should return a vector of shape (n_samples,) of type int or float, \n default=None\n function used to compute Mondrian categories\n de : a :class:`.DifficultyEstimator`, default=None\n a fitted difficulty estimator (used only if f is not None)\n learner : an object with the method ``learner.predict``, default=None\n a fitted regression model (used only if de and f are not None) \n no_bins : int, default=10\n no. of Mondrian categories\n oob : bool, default=False\n use out-of-bag estimation (not used if f is not None)\n\n Returns\n -------\n self : object\n Fitted MondrianCategorizer.\n\n Examples\n --------\n Assuming that ``X_train`` is an array of shape (n_samples, n_features)\n and ``get_values`` is a function that given ``X_train`` returns a vector\n of values of shape (n_samples,), then a Mondrian categorizer can\n be formed in the following way, where the boundaries for the Mondrian\n categories are found by partitioning the values in the vector into five\n equal-sized bins:\n \n .. code-block:: python\n\n from crepes.extras import MondrianCategorizer\n\n mc = MondrianCategorizer()\n mc.fit(X, f=get_values, no_bins=5)\n \"\"\"\n if f is not None:\n if X is not None:\n scores = f(X)\n bins, bin_thresholds = binning(scores, bins=no_bins)\n self.bin_thresholds = bin_thresholds\n else:\n raise ValueError(\"X must be provided since f is not None\")\n self.f = f\n elif de is not None:\n if oob:\n scores = de.apply()\n self.oob = True\n else:\n if X is not None:\n scores = de.apply(X)\n else:\n raise ValueError((\"X must be provided since de is not None\"\n \"and oob=False\"))\n self.de = de\n bins, bin_thresholds = binning(scores, bins=no_bins)\n self.bin_thresholds = bin_thresholds \n elif learner is not None:\n if oob:\n scores = learner.oob_prediction_\n self.oob = True\n else:\n if X is not None:\n scores = learner.predict(X)\n else:\n raise ValueError((\"X must be provided since learner is not None\"\n \"and oob=False\"))\n self.learner = learner\n bins, bin_thresholds = binning(scores, bins=no_bins)\n self.bin_thresholds = bin_thresholds\n else:\n raise ValueError(\"One of f, de, and learner must not be None\")\n self.fitted = True\n return self\n\n def apply(self, X):\n \"\"\"\n Apply Mondrian categorizer.\n\n Parameters\n ----------\n X : array-like of shape (n_samples, n_features)\n set of objects\n\n Returns\n -------\n bins : array-like of shape (n_samples,)\n Mondrian categories \n\n Examples\n --------\n Assuming ``mc`` to be a fitted :class:`.MondrianCategorizer`, i.e., for\n which :meth:`.fit` has earlier been called, Mondrian categories for a\n set of objects ``X`` is obtained by:\n\n .. code-block:: python\n \n categories = mc.apply(X)\n\n Note\n ----\n The array used when calling :meth:`.fit` must have the same number of\n columns (``n_features``) as the array used as input to :meth:`.apply`.\n \"\"\"\n if self.f is not None:\n if self.bin_thresholds is None:\n bins = self.f(X)\n else:\n scores = self.f(X)\n bins = binning(scores, bins=self.bin_thresholds)\n elif self.de is not None:\n scores = self.de.apply(X)\n bins = binning(scores, bins=self.bin_thresholds)\n elif self.learner is not None:\n if self.oob:\n predictions = np.array([model.predict(X)\n for model in self.learner.estimators_])\n oob_masks = np.array([\n get_oob(self.learner.estimators_[i].random_state, len(X))\n for i in range(len(self.learner.estimators_))])\n scores = np.array([np.mean(predictions[oob_masks[:,i],i])\n for i in range(len(X))])\n else:\n scores = learner.predict(X)\n bins = binning(scores, bins=self.bin_thresholds)\n return bins\n \nclass DifficultyEstimator():\n \"\"\"\n A difficulty estimator outputs scores for objects to be used by \n normalized conformal regressors and predictive systems.\n \"\"\"\n \n def __init__(self):\n self.fitted = False\n self.estimator_type = None\n\n def __repr__(self):\n if self.fitted and self.estimator_type == \"knn\":\n return (f\"DifficultyEstimator(fitted={self.fitted}, \"\n f\"type={self.estimator_type}, \"\n f\"k={self.k}, \"\n f\"target={self.target_type}, \"\n f\"scaler={self.scaler}, \" \n f\"beta={self.beta}, \"\n f\"oob={self.oob}\" \n \")\")\n elif self.fitted and self.estimator_type == \"variance\":\n return (f\"DifficultyEstimator(fitted={self.fitted}, \"\n f\"type={self.estimator_type}, \"\n f\"scaler={self.scaler}, \" \n f\"beta={self.beta}, \"\n f\"oob={self.oob}\" \n \")\")\n else:\n return f\"DifficultyEstimator(fitted={self.fitted})\"\n \n def fit(self, X=None, f=None, y=None, residuals=None, learner=None,\n k=25, scaler=False, beta=0.01, oob=False):\n \"\"\"\n Fit difficulty estimator.\n\n Parameters\n ----------\n X : array-like of shape (n_samples, n_features), default=None\n set of objects\n f : function which given an array-like of shape (n_samples, n_features)\n should return a vector of shape (n_samples,) of type int or float, \n default=None\n function used to compute difficulty estimates\n y : array-like of shape (n_samples,), default=None\n target values\n residuals : array-like of shape (n_samples,), default=None\n true target values - predicted values\n learner : an object with attribute ``learner.estimators_``, default=None\n an ensemble model where each model m in ``learner.estimators_`` has a\n method ``m.predict`` (used only if f=None)\n k: int, default=25\n number of neighbors (used only if f=None and learner=None)\n scaler : bool, default=True\n use min-max-scaler on the difficulty estimates\n beta : int or float, default=0.01 \n value to add to the difficulty estimates (after scaling)\n oob : bool, default=False\n use out-of-bag estimation\n\n Returns\n -------\n self : object\n Fitted DifficultyEstimator.\n\n Examples\n --------\n Assuming that ``X_prop_train`` is a proper training set, \n then a difficulty estimator using the distances to the k \n nearest neighbors can be formed in the following way \n (here using the default ``k=25``):\n \n .. code-block:: python\n\n from crepes.extras import DifficultyEstimator\n\n de_knn_dist = DifficultyEstimator() \n de_knn_dist.fit(X_prop_train)\n\n Assuming that ``y_prop_train`` is a vector with target values \n for the proper training set, then a difficulty estimator using \n standard deviation of the targets of the k nearest neighbors \n is formed by: \n\n .. code-block:: python\n\n de_knn_std = DifficultyEstimator() \n de_knn_std.fit(X_prop_train, y=y_prop_train)\n\n Assuming that ``X_prop_res`` is a vector with residuals \n for the proper training set, then a difficulty estimator using \n the mean of the absolute residuals of the k nearest neighbors \n is formed by: \n\n .. code-block:: python\n\n de_knn_res = DifficultyEstimator() \n de_knn_res.fit(X_prop_train, residuals=X_prop_res)\n\n Assuming that ``learner_prop`` is a trained model for which\n ``learner.estimators_`` is a collection of base models, each \n implementing the ``predict`` method; this holds e.g., for \n ``RandomForestRegressor``, a difficulty estimator using the variance\n of the predictions of the constituent models is formed by: \n\n .. code-block:: python\n\n de_var = DifficultyEstimator() \n de_var.fit(learner=learner_prop)\n\n The difficulty estimates may be normalized (using min-max scaling) by\n setting ``scaler=True``. It should be noted that this comes with a \n computational cost; for estimators based on the k-nearest neighbor, \n a leave-one-out protocol is employed to find the minimum and maximum \n distances that are used by the scaler. This also requires that a set \n of objects is provided for the variance-based approach (to allow for \n finding the minimum and maximum values). Hence, if normalization is to\n be employed for the latter, objects have to be included:\n\n .. code-block:: python\n\n de_var = DifficultyEstimator() \n de_var.fit(X_proper_train, learner=learner_prop, scaler=True)\n\n Difficulty estimates may also be computed by an externally defined\n function. Assuming that ``diff_model`` is a fitted regression model,\n for which the ``predict`` method gives estimates of the absolute\n error for the objects in ``X_proper_train``, then normalized difficulty\n estimates can be obtained from the following difficulty estimator:\n\n .. code-block:: python\n\n de_mod = DifficultyEstimator() \n de_mod.fit(X_proper_train, f=diff_model.predict, scaler=True)\n \n The :class:`.DifficultyEstimator` can also support the construction of \n conformal regressors and predictive systems that employ out-of-bag \n calibration. For the k-nearest neighbor approaches, the difficulty of\n each object in the provided training set will be computed using a \n leave-one-out procedure, while for the variance-based approach the \n out-of-bag predictions will be employed. This is enabled by setting \n ``oob=True`` when calling the :meth:`.fit` method, which also requires \n the (full) training set (``X_train``), and for the variance-based \n approach a corresponding trained model (``learner_full``) to be \n provided: \n\n .. code-block:: python\n\n de_var_oob = DifficultyEstimator() \n de_var_oob.fit(X_train, learner=learner_full, scaler=True, oob=True)\n\n A small value (beta) is added to the difficulty estimates. The default \n is ``beta=0.01``. In order to make the beta value have the same effect \n across different estimators, you may consider normalizing the difficulty\n estimates (using min-max scaling) by setting ``scaler=True``. Note that \n beta is added after the normalization, which means that the range of\n scores after normalization will be [0+beta, 1+beta]. Below, we use \n ``beta=0.001`` together with 10 neighbors (``k=10``):\n\n .. code-block:: python\n \n de_knn_mod = DifficultyEstimator() \n de_knn_mod.fit(X_prop_train, k=10, beta=0.001, scaler=True)\n\n Note\n ----\n The use of out-of-bag calibration, as enabled by ``oob=True``, \n does not come with the theoretical validity guarantees of the regular\n (inductive) conformal regressors and predictive systems, due to that \n calibration and test instances are not handled in exactly the same way.\n \"\"\"\n self.f = f\n if isinstance(y, pd.Series):\n y = y.values\n self.y = y\n self.residuals = residuals\n self.learner = learner\n self.k = k\n self.beta = beta\n self.scaler = scaler\n self.oob = oob\n if self.f is not None:\n self.estimator_type = \"function\"\n elif self.learner is not None:\n self.estimator_type = \"variance\"\n try:\n self.learner.estimators_\n except:\n raise ValueError(\n \"learner is missing the attribute estimators_\")\n if self.oob:\n try:\n self.learner.estimators_[0].random_state\n except:\n raise ValueError(\n (\"learner.estimators_ is missing the attribute \"\n \"random_state\"))\n else:\n self.estimator_type = \"knn\"\n if self.residuals is None:\n if self.y is None:\n self.target_type = \"none\"\n else:\n self.target_type = \"labels\"\n else:\n self.target_type = \"residuals\"\n if X is None:\n raise ValueError(\"X=None is not allowed for k-nearest\"\n \" neighbor estimators\")\n \n if self.estimator_type == \"function\":\n if X is None and self.scaler:\n raise ValueError(\"X=None is allowed only if scaler=False\"\n \" for function estimators\")\n if self.oob:\n raise ValueError(\"oob=True is not allowed for function\"\n \" estimators\")\n if self.scaler:\n sigmas = self.f(X)\n sigma_scaler = MinMaxScaler(clip=True)\n sigma_scaler.fit(sigmas[:,None])\n self.sigma_scaler = sigma_scaler\n\n if self.estimator_type == \"variance\":\n if X is None and (self.oob or self.scaler):\n raise ValueError(\"X=None is allowed only if oob=False and \"\n \"scaler=False for variance estimator\")\n if self.oob or self.scaler:\n predictions = np.array([model.predict(X)\n for model in self.learner.estimators_])\n oob_masks = np.array([\n get_oob(self.learner.estimators_[i].random_state, len(X))\n for i in range(len(self.learner.estimators_))])\n sigmas = np.array([np.var(predictions[oob_masks[:,i],i])\n for i in range(len(X))])\n if self.scaler:\n sigma_scaler = MinMaxScaler(clip=True)\n sigma_scaler.fit(sigmas[:,None])\n self.sigma_scaler = sigma_scaler\n sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] \n if self.oob:\n self.sigmas = sigmas\n\n if self.estimator_type == \"knn\":\n nn = NearestNeighbors(n_neighbors=self.k, n_jobs=-1)\n nn_scaler = MinMaxScaler(clip=True)\n nn_scaler.fit(X)\n X_scaled = nn_scaler.transform(X)\n nn.fit(X_scaled)\n self.nn = nn\n self.nn_scaler = nn_scaler\n if self.oob or self.scaler:\n if self.target_type == \"none\":\n distances, neighbor_indexes = nn.kneighbors(\n return_distance=True)\n sigmas = np.array([np.sum(distances[i])\n for i in range(len(distances))]) \n elif self.target_type == \"labels\":\n neighbor_indexes = nn.kneighbors(return_distance=False)\n sigmas = np.array([np.std(y[indexes])\n for indexes in neighbor_indexes])\n else: # self.target_type == \"residuals\"\n neighbor_indexes = nn.kneighbors(return_distance=False)\n sigmas = np.array([np.mean(np.abs(residuals[indexes]))\n for indexes in neighbor_indexes])\n if self.scaler:\n sigma_scaler = MinMaxScaler(clip=True)\n sigma_scaler.fit(sigmas[:,None])\n self.sigma_scaler = sigma_scaler\n sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0]\n if self.oob:\n self.sigmas = sigmas\n \n self.fitted = True\n return self\n\n def apply(self, X=None):\n \"\"\"\n Apply difficulty estimator.\n\n Parameters\n ----------\n X : array-like of shape (n_samples, n_features), default=None\n set of objects\n\n Returns\n -------\n sigmas : array-like of shape (n_samples,)\n difficulty estimates \n\n Examples\n --------\n Assuming ``de`` to be a fitted :class:`.DifficultyEstimator`, i.e., for\n which :meth:`.fit` has earlier been called, difficulty estimates for a\n set of objects ``X`` is obtained by:\n\n .. code-block:: python\n \n difficulty_estimates = de.apply(X)\n\n If ``de_oob`` is a :class:`.DifficultyEstimator` that has been fitted \n with the option ``oob=True`` and a training set, then a call to \n :meth:`.apply` without any objects will return the estimates for the \n training set:\n\n .. code-block:: python\n \n oob_difficulty_estimates = de_oob.apply()\n\n For a difficulty estimator employing any of the k-nearest neighbor \n approaches, the above will return an estimate for the difficulty \n of each object in the training set computed using a leave-one-out \n procedure, while for the variance-based approach the out-of-bag \n predictions will instead be used. \n \"\"\"\n if X is None:\n if not self.oob:\n raise ValueError(\"X=None is allowed only if oob=True\")\n sigmas = self.sigmas\n elif self.estimator_type == \"knn\":\n X_scaled = self.nn_scaler.transform(X)\n if self.target_type == \"none\":\n distances, neighbor_indexes = self.nn.kneighbors(\n X_scaled, return_distance=True)\n sigmas = np.array([np.sum(distances[i])\n for i in range(len(distances))])\n elif self.target_type == \"labels\":\n neighbor_indexes = self.nn.kneighbors(X_scaled,\n return_distance=False)\n sigmas = np.array([np.std(self.y[indexes])\n for indexes in neighbor_indexes])\n else: # self.target_type == \"residuals\"\n neighbor_indexes = self.nn.kneighbors(X_scaled,\n return_distance=False)\n sigmas = np.array([np.mean(np.abs(self.residuals[indexes]))\n for indexes in neighbor_indexes])\n if self.scaler:\n sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0]\n elif self.estimator_type == \"variance\":\n if self.oob:\n predictions = np.array([model.predict(X)\n for model in self.learner.estimators_])\n oob_masks = np.array([\n get_oob(self.learner.estimators_[i].random_state, len(X))\n for i in range(len(self.learner.estimators_))])\n sigmas = np.array([np.var(predictions[oob_masks[:,i],i])\n for i in range(len(X))])\n else: \n sigmas = np.var([model.predict(X) for\n model in self.learner.estimators_],\n axis=0)\n if self.scaler:\n sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] \n else: # self.estimator_type == \"function\"\n sigmas = self.f(X)\n if self.scaler:\n sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] \n return sigmas + self.beta\n\ndef get_oob(seed, n_samples):\n \"\"\"\n Provides out-of-bag samples from a random seed and sample size.\n\n Parameters\n ----------\n seed : int\n random seed\n n_samples : int\n sample size\n \n Returns\n -------\n oob : array-like of shape (n_samples,)\n binary vector indicating which samples are out-of-bag and not \n \"\"\"\n return np.bincount(np.random.RandomState(seed).randint(0, n_samples,\n n_samples),\n minlength=n_samples) == 0","metadata":{"trusted":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2024-11-06T01:51:26.334579Z","iopub.execute_input":"2024-11-06T01:51:26.335171Z","iopub.status.idle":"2024-11-06T01:51:26.430310Z","shell.execute_reply.started":"2024-11-06T01:51:26.335103Z","shell.execute_reply":"2024-11-06T01:51:26.429080Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"Conformal classifiers, regressors, and predictive systems (crepes)\n\nClasses implementing conformal classifiers, regressors, and predictive\nsystems, on top of any standard classifier and regressor, transforming\nthe original predictions into well-calibrated p-values and cumulative\ndistribution functions, or prediction sets and intervals with coverage\nguarantees.\n\nAuthor: Henrik Boström (bostromh@kth.se)\n\nCopyright 2024 Henrik Boström\n\nLicense: BSD 3 clause\n\n\"\"\"\n\nimport numpy as np\nimport pandas as pd\nimport time\n\nclass ConformalPredictor():\n \"\"\"\n The class contains three sub-classes: :class:`.ConformalClassifier`,\n :class:`.ConformalRegressor`, and :class:`.ConformalPredictiveSystem`.\n \"\"\"\n \n def __init__(self):\n self.fitted = False\n self.mondrian = None\n self.time_fit = None\n self.time_predict = None\n self.time_evaluate = None\n self.alphas = None\n self.normalized = None\n self.seed = None\n\nclass ConformalClassifier(ConformalPredictor):\n \"\"\"\n A conformal classifier transforms non-conformity scores into p-values\n or prediction sets for a certain confidence level.\n \"\"\"\n \n def __repr__(self):\n if self.fitted:\n return (f\"ConformalClassifier(fitted={self.fitted}, \"\n f\"mondrian={self.mondrian})\")\n else:\n return f\"ConformalClassifier(fitted={self.fitted})\"\n \n def fit(self, alphas, bins=None, seed=None):\n \"\"\"\n Fit conformal classifier.\n\n Parameters\n ----------\n alphas : array-like of shape (n_samples,)\n non-conformity scores\n bins : array-like of shape (n_samples,), default=None\n Mondrian categories\n seed : int, default=None\n set random seed\n\n Returns\n -------\n self : object\n Fitted ConformalClassifier.\n\n Examples\n --------\n Assuming that ``alphas_cal`` is a vector with non-conformity scores,\n then a standard conformal classifier is formed in the following way:\n\n .. code-block:: python\n\n from crepes import ConformalClassifier\n\n cc_std = ConformalClassifier() \n\n cc_std.fit(alphas_cal) \n\n Assuming that ``bins_cals`` is a vector with Mondrian categories \n (bin labels), then a Mondrian conformal classifier is fitted in the\n following way:\n\n .. code-block:: python\n\n cc_mond = ConformalClassifier()\n cc_mond.fit(alphas_cal, bins=bins_cal)\n\n Note\n ----\n By providing a random seed, e.g., ``seed=123``, calls to the methods\n ``predict_p``, ``predict_set`` and ``evaluate`` of the\n :class:`.ConformalClassifier` object will be deterministic.\n \"\"\"\n tic = time.time()\n if bins is None:\n self.mondrian = False\n self.alphas = np.sort(alphas)[::-1]\n else: \n self.mondrian = True\n bin_values = np.unique(bins)\n self.alphas = (bin_values, [np.sort(alphas[bins==b])[::-1]\n for b in bin_values])\n self.seed = seed\n self.fitted = True\n toc = time.time()\n self.time_fit = toc-tic\n return self\n\n def predict_p(self, alphas, bins=None, confidence=0.95, smoothing=True,\n seed=None):\n \"\"\"\n Obtain (smoothed or non-smoothed) p-values from conformal classifier.\n\n Parameters\n ----------\n alphas : array-like of shape (n_samples, n_classes)\n non-conformity scores\n bins : array-like of shape (n_samples,), default=None\n Mondrian categories\n confidence : float in range (0,1), default=0.95\n confidence level\n smoothing : bool, default=True\n use smoothed p-values\n seed : int, default=None\n set random seed\n\n Returns\n -------\n p-values : ndarray of shape (n_samples, n_classes)\n p-values\n\n Examples\n --------\n Assuming that ``alphas_test`` is a vector with non-conformity scores\n for a test set and ``cc_std`` a fitted standard conformal classifier, \n then p-values for the test is obtained by:\n\n .. code-block:: python\n\n p_values = cc_std.predict_p(alphas_test)\n\n Assuming that ``bins_test`` is a vector with Mondrian categories (bin \n labels) for the test set and ``cc_mond`` a fitted Mondrian conformal \n classifier, then the following provides (smoothed) p-values for the\n test set:\n\n .. code-block:: python\n\n p_values = cc_mond.predict_p(alphas_test, bins=bins_test)\n\n Note\n ----\n If a value for ``seed`` is given, it will take precedence over any ``seed``\n value given when calling ``fit``.\n \"\"\"\n tic = time.time()\n if seed is None:\n seed = self.seed\n if seed is not None:\n random_state = np.random.get_state()\n np.random.seed(seed)\n if not self.mondrian:\n if smoothing:\n p_values = np.array(\n [[(np.sum(self.alphas > alpha) + np.random.rand()*(\n np.sum(self.alphas == alpha)+1))/(len(self.alphas)+1)\n for alpha in alpha_row] for alpha_row in alphas])\n else:\n p_values = np.array(\n [[(np.sum(self.alphas >= alpha) + 1)/(len(self.alphas)+1)\n for alpha in alpha_row] for alpha_row in alphas])\n else:\n bin_values, bin_alphas = self.alphas\n bin_indexes = np.array([np.argwhere(bin_values == bins[i])[0][0]\n for i in range(len(bins))])\n if smoothing:\n p_values = np.array([\n [(np.sum(bin_alphas[bin_indexes[i]] > alpha) \\\n + np.random.rand()*(np.sum(bin_alphas[\n bin_indexes[i]] == alpha)+1))/(\n len(bin_alphas[bin_indexes[i]])+1)\n for alpha in alphas[i]]\n for i in range(len(alphas))])\n else:\n p_values = np.array([\n [(np.sum(bin_alphas[bin_indexes[i]] >= alpha) + 1)/(\n len(bin_alphas[bin_indexes[i]])+1)\n for alpha in alphas[i]]\n for i in range(len(alphas))])\n if seed is not None:\n np.random.set_state(random_state)\n toc = time.time()\n self.time_predict = toc-tic \n return p_values\n \n def predict_set(self, alphas, bins=None, confidence=0.95, smoothing=True,\n seed=None):\n \"\"\"\n Obtain prediction sets using conformal classifier.\n\n Parameters\n ----------\n alphas : array-like of shape (n_samples, n_classes)\n non-conformity scores\n bins : array-like of shape (n_samples,), default=None\n Mondrian categories\n confidence : float in range (0,1), default=0.95\n confidence level\n smoothing : bool, default=True\n use smoothed p-values\n seed : int, default=None\n set random seed\n\n Returns\n -------\n prediction sets : ndarray of shape (n_samples, n_classes)\n prediction sets\n\n Examples\n --------\n Assuming that ``alphas_test`` is a vector with non-conformity scores\n for a test set and ``cc_std`` a fitted standard conformal classifier, \n then prediction sets at the default (95%) confidence level are\n obtained by:\n\n .. code-block:: python\n\n prediction_sets = cc_std.predict_set(alphas_test)\n\n Assuming that ``bins_test`` is a vector with Mondrian categories (bin \n labels) for the test set and ``cc_mond`` a fitted Mondrian conformal \n classifier, then the following provides prediction sets for the test set,\n at the 90% confidence level:\n\n .. code-block:: python\n\n p_values = cc_mond.predict_set(alphas_test, \n bins=bins_test,\n confidence=0.9)\n\n Note\n ----\n The use of smoothed p-values increases computation time and typically\n has a minor effect on the predictions sets, except for small calibration\n sets.\n\n Note\n ----\n If a value for ``seed`` is given, it will take precedence over any ``seed``\n value given when calling ``fit``.\n \"\"\"\n tic = time.time()\n if seed is None:\n seed = self.seed\n if seed is not None:\n random_state = np.random.get_state()\n np.random.seed(seed)\n if smoothing:\n p_values = self.predict_p(alphas, bins, smoothing=True)\n prediction_sets = (p_values >= 1-confidence).astype(int)\n elif bins is None:\n alpha_index = int((1-confidence)*(len(self.alphas)+1))-1\n if alpha_index >= 0:\n alpha_value = self.alphas[alpha_index]\n prediction_sets = (alphas <= alpha_value).astype(int)\n else:\n prediction_sets = np.ones(alphas.shape)\n warnings.warn(\"the no. of calibration examples is \" \\\n \"too small for the chosen confidence level; \" \\\n \"all labels are included in the prediction sets\")\n else:\n bin_values, bin_alphas = self.alphas\n alpha_indexes = np.array(\n [int((1-confidence)*(len(bin_alphas[b])+1))-1\n for b in range(len(bin_values))])\n alpha_values = [bin_alphas[b][alpha_indexes[b]]\n if alpha_indexes[b] >= 0\n else -np.inf for b in range(len(bin_values))]\n bin_indexes = np.array([np.argwhere(bin_values == bins[i])[0][0]\n for i in range(len(bins))])\n prediction_sets = np.array(\n [alphas[i] <= alpha_values[bin_indexes[i]]\n for i in range(len(alphas))], dtype=int)\n if (alpha_indexes < 0).any():\n warnings.warn(\"the no. of calibration examples in some bins is\" \\\n \" too small for the chosen confidence level; \" \\\n \"all labels are included in the corresponding\" \\\n \"prediction sets\")\n if seed is not None:\n np.random.set_state(random_state)\n toc = time.time()\n self.time_predict = toc-tic \n return prediction_sets\n\n def evaluate(self, alphas, classes, y, bins=None, confidence=0.95,\n smoothing=True, metrics=None, seed=None):\n \"\"\"\n Evaluate conformal classifier.\n\n Parameters\n ----------\n alphas : array-like of shape (n_samples, n_classes)\n non-conformity scores\n classes : array-like of shape (n_classes,)\n class names\n y : array-like of shape (n_samples,)\n correct class labels\n bins : array-like of shape (n_samples,), default=None\n Mondrian categories\n confidence : float in range (0,1), default=0.95\n confidence level\n smoothing : bool, default=True\n use smoothed p-values\n metrics : a string or a list of strings, \n default = list of all metrics, i.e., [\"error\", \"avg_c\", \n \"one_c\", \"empty\", \"time_fit\", \"time_evaluate\"]\n seed : int, default=None\n set random seed\n \n Returns\n -------\n results : dictionary with a key for each selected metric \n estimated performance using the metrics, where \"error\" is the \n fraction of prediction sets not containing the true class label,\n \"avg_c\" is the average no. of predicted class labels, \"one_c\" is\n the fraction of singleton prediction sets, \"empty\" is the fraction\n of empty prediction sets, \"time_fit\" is the time taken to fit the \n conformal classifier, and \"time_evaluate\" is the time taken for the\n evaluation \n\n Examples\n --------\n Assuming that ``alphas`` is an array containing non-conformity scores \n for all classes for the test objects, ``classes`` and ``y_test`` are \n vectors with the class names and true class labels for the test set, \n respectively, and ``cc`` is a fitted standard conformal classifier, \n then the latter can be evaluated at the default confidence level with \n respect to error and average number of labels in the prediction sets by:\n\n .. code-block:: python\n\n results = cc.evaluate(alphas, y_test, metrics=[\"error\", \"avg_c\"])\n\n Note\n ----\n The use of smoothed p-values increases computation time and typically\n has a minor effect on the results, except for small calibration sets.\n\n Note\n ----\n If a value for ``seed`` is given, it will take precedence over any ``seed``\n value given when calling ``fit``. \n \"\"\"\n if metrics is None:\n metrics = [\"error\", \"avg_c\", \"one_c\", \"empty\", \"time_fit\",\n \"time_evaluate\"]\n tic = time.time()\n if seed is None:\n seed = self.seed\n if seed is not None:\n random_state = np.random.get_state()\n np.random.seed(seed)\n prediction_sets = self.predict_set(alphas, bins, confidence, smoothing)\n test_results = get_test_results(prediction_sets, classes, y, metrics)\n if seed is not None:\n np.random.set_state(random_state)\n toc = time.time()\n self.time_evaluate = toc-tic\n if \"time_fit\" in metrics:\n test_results[\"time_fit\"] = self.time_fit\n if \"time_evaluate\" in metrics:\n test_results[\"time_evaluate\"] = self.time_evaluate\n return test_results\n\ndef get_test_results(prediction_sets, classes, y, metrics):\n test_results = {}\n class_indexes = np.array(\n [np.argwhere(classes == y[i])[0][0] for i in range(len(y))]) \n if \"error\" in metrics:\n test_results[\"error\"] = 1-np.sum(\n prediction_sets[np.arange(len(y)), class_indexes]) / len(y)\n if \"avg_c\" in metrics: \n test_results[\"avg_c\"] = np.sum(prediction_sets) / len(y)\n if \"one_c\" in metrics: \n test_results[\"one_c\"] = np.sum(\n [np.sum(p) == 1 for p in prediction_sets]) / len(y)\n if \"empty\" in metrics: \n test_results[\"empty\"] = np.sum(\n [np.sum(p) == 0 for p in prediction_sets]) / len(y)\n return test_results\n\nclass ConformalRegressor(ConformalPredictor):\n \"\"\"\n A conformal regressor transforms point predictions (regression \n values) into prediction intervals, for a certain confidence level.\n \"\"\"\n \n def __repr__(self):\n if self.fitted:\n return (f\"ConformalRegressor(fitted={self.fitted}, \"\n f\"normalized={self.normalized}, \"\n f\"mondrian={self.mondrian})\")\n else:\n return f\"ConformalRegressor(fitted={self.fitted})\"\n \n def fit(self, residuals, sigmas=None, bins=None):\n \"\"\"\n Fit conformal regressor.\n\n Parameters\n ----------\n residuals : array-like of shape (n_values,)\n true values - predicted values\n sigmas: array-like of shape (n_values,), default=None\n difficulty estimates\n bins : array-like of shape (n_values,), default=None\n Mondrian categories\n\n Returns\n -------\n self : object\n Fitted ConformalRegressor.\n\n Examples\n --------\n Assuming that ``y_cal`` and ``y_hat_cal`` are vectors with true\n and predicted targets for some calibration set, then a standard\n conformal regressor can be formed from the residuals:\n\n .. code-block:: python\n\n residuals_cal = y_cal - y_hat_cal\n\n from crepes import ConformalRegressor\n\n cr_std = ConformalRegressor() \n\n cr_std.fit(residuals_cal) \n\n Assuming that ``sigmas_cal`` is a vector with difficulty estimates,\n then a normalized conformal regressor can be fitted in the following\n way:\n\n .. code-block:: python\n\n cr_norm = ConformalRegressor()\n cr_norm.fit(residuals_cal, sigmas=sigmas_cal)\n\n Assuming that ``bins_cals`` is a vector with Mondrian categories \n (bin labels), then a Mondrian conformal regressor can be fitted in the\n following way:\n\n .. code-block:: python\n\n cr_mond = ConformalRegressor()\n cr_mond.fit(residuals_cal, bins=bins_cal)\n\n A normalized Mondrian conformal regressor can be fitted in the \n following way:\n\n .. code-block:: python\n\n cr_norm_mond = ConformalRegressor()\n cr_norm_mond.fit(residuals_cal, sigmas=sigmas_cal, \n bins=bins_cal)\n \"\"\"\n tic = time.time()\n abs_residuals = np.abs(residuals)\n if bins is None:\n self.mondrian = False\n if sigmas is None:\n self.normalized = False\n self.alphas = np.sort(abs_residuals)[::-1]\n else:\n self.normalized = True\n self.alphas = np.sort(abs_residuals/sigmas)[::-1]\n else: \n self.mondrian = True\n bin_values = np.unique(bins)\n if sigmas is None: \n self.normalized = False\n self.alphas = (bin_values,[np.sort(\n abs_residuals[bins==b])[::-1] for b in bin_values])\n else:\n self.normalized = True\n self.alphas = (bin_values, [np.sort(\n abs_residuals[bins==b]/sigmas[bins==b])[::-1]\n for b in bin_values]) \n self.fitted = True\n toc = time.time()\n self.time_fit = toc-tic\n return self\n\n def predict(self, y_hat, sigmas=None, bins=None, confidence=0.95,\n y_min=-np.inf, y_max=np.inf):\n \"\"\"\n Predict using conformal regressor.\n\n Parameters\n ----------\n y_hat : array-like of shape (n_values,)\n predicted values\n sigmas : array-like of shape (n_values,), default=None\n difficulty estimates\n bins : array-like of shape (n_values,), default=None\n Mondrian categories\n confidence : float in range (0,1), default=0.95\n confidence level\n y_min : float or int, default=-numpy.inf\n minimum value to include in prediction intervals\n y_max : float or int, default=numpy.inf\n maximum value to include in prediction intervals\n\n Returns\n -------\n intervals : ndarray of shape (n_values, 2)\n prediction intervals\n\n Examples\n --------\n Assuming that ``y_hat_test`` is a vector with predicted targets for a\n test set and ``cr_std`` a fitted standard conformal regressor, then \n prediction intervals at the 99% confidence level can be obtained by:\n\n .. code-block:: python\n\n intervals = cr_std.predict(y_hat_test, confidence=0.99)\n\n Assuming that ``sigmas_test`` is a vector with difficulty estimates for\n the test set and ``cr_norm`` a fitted normalized conformal regressor, \n then prediction intervals at the default (95%) confidence level can be\n obtained by:\n\n .. code-block:: python\n\n intervals = cr_norm.predict(y_hat_test, sigmas=sigmas_test)\n\n Assuming that ``bins_test`` is a vector with Mondrian categories (bin \n labels) for the test set and ``cr_mond`` a fitted Mondrian conformal \n regressor, then the following provides prediction intervals at the \n default confidence level, where the intervals are lower-bounded by 0:\n\n .. code-block:: python\n\n intervals = cr_mond.predict(y_hat_test, bins=bins_test, \n y_min=0)\n\n Note\n ----\n In case the specified confidence level is too high in relation to the \n size of the calibration set, a warning will be issued and the output\n intervals will be of maximum size.\n \"\"\"\n tic = time.time()\n intervals = np.zeros((len(y_hat),2))\n if not self.mondrian:\n alpha_index = int((1-confidence)*(len(self.alphas)+1))-1\n if alpha_index >= 0:\n alpha = self.alphas[alpha_index]\n if self.normalized:\n intervals[:,0] = y_hat-alpha*sigmas\n intervals[:,1] = y_hat+alpha*sigmas\n else:\n intervals[:,0] = y_hat-alpha\n intervals[:,1] = y_hat+alpha\n else:\n intervals[:,0] = -np.inf \n intervals[:,1] = np.inf\n warnings.warn(\"the no. of calibration examples is too small\" \\\n \"for the chosen confidence level; the \" \\\n \"intervals will be of maximum size\")\n else: \n bin_values, bin_alphas = self.alphas\n bin_indexes = [np.argwhere(bins == b).T[0]\n for b in bin_values]\n alpha_indexes = np.array(\n [int((1-confidence)*(len(bin_alphas[b])+1))-1\n for b in range(len(bin_values))])\n too_small_bins = np.argwhere(alpha_indexes < 0)\n if len(too_small_bins) > 0:\n if len(too_small_bins[:,0]) < 11:\n bins_to_show = \" \".join([str(bin_values[i]) for i in\n too_small_bins[:,0]])\n else:\n bins_to_show = \" \".join([str(bin_values[i]) for i in\n too_small_bins[:10,0]]+['...'])\n warnings.warn(\"the no. of calibration examples is too \" \\\n \"small for the chosen confidence level \" \\\n f\"in the following bins: {bins_to_show}; \"\\\n \"the corresponding intervals will be of \" \\\n \"maximum size\") \n bin_alpha = np.array([bin_alphas[b][alpha_indexes[b]]\n if alpha_indexes[b]>=0 else np.inf\n for b in range(len(bin_values))])\n if self.normalized:\n for b in range(len(bin_values)):\n intervals[bin_indexes[b],0] = y_hat[bin_indexes[b]] \\\n - bin_alpha[b]*sigmas[bin_indexes[b]]\n intervals[bin_indexes[b],1] = y_hat[bin_indexes[b]] \\\n + bin_alpha[b]*sigmas[bin_indexes[b]]\n else:\n for b in range(len(bin_values)):\n intervals[bin_indexes[b],0] = y_hat[bin_indexes[b]] \\\n - bin_alpha[b]\n intervals[bin_indexes[b],1] = y_hat[bin_indexes[b]] \\\n + bin_alpha[b] \n if y_min > -np.inf:\n intervals[intervalsy_max] = y_max \n toc = time.time()\n self.time_predict = toc-tic \n return intervals\n\n def evaluate(self, y_hat, y, sigmas=None, bins=None,\n confidence=0.95, y_min=-np.inf, y_max=np.inf, metrics=None):\n \"\"\"\n Evaluate conformal regressor.\n\n Parameters\n ----------\n y_hat : array-like of shape (n_values,)\n predicted values\n y : array-like of shape (n_values,)\n correct target values\n sigmas : array-like of shape (n_values,), default=None\n difficulty estimates\n bins : array-like of shape (n_values,), default=None\n Mondrian categories\n confidence : float in range (0,1), default=0.95\n confidence level\n y_min : float or int, default=-numpy.inf\n minimum value to include in prediction intervals\n y_max : float or int, default=numpy.inf\n maximum value to include in prediction intervals\n metrics : a string or a list of strings, \n default=list of all metrics, i.e., \n [\"error\", \"eff_mean\", \"eff_med\", \"time_fit\", \"time_evaluate\"]\n \n Returns\n -------\n results : dictionary with a key for each selected metric \n estimated performance using the metrics\n\n Examples\n --------\n Assuming that ``y_hat_test`` and ``y_test`` are vectors with predicted\n and true targets for a test set, ``sigmas_test`` and ``bins_test`` are\n vectors with difficulty estimates and Mondrian categories (bin labels) \n for the test set, and ``cr_norm_mond`` is a fitted normalized Mondrian\n conformal regressor, then the latter can be evaluated at the default\n confidence level with respect to error and mean efficiency (interval \n size) by:\n\n .. code-block:: python\n\n results = cr_norm_mond.evaluate(y_hat_test, y_test, \n sigmas=sigmas_test, bins=bins_test,\n metrics=[\"error\", \"eff_mean\"])\n \"\"\"\n tic = time.time()\n if metrics is None:\n metrics = [\"error\",\"eff_mean\",\"eff_med\",\"time_fit\",\"time_evaluate\"]\n test_results = {}\n intervals = self.predict(y_hat, sigmas, bins, confidence, y_min, y_max)\n if \"error\" in metrics:\n test_results[\"error\"] = 1-np.mean(\n np.logical_and(intervals[:,0]<=y, y<=intervals[:,1]))\n if \"eff_mean\" in metrics: \n test_results[\"eff_mean\"] = np.mean(intervals[:,1]-intervals[:,0])\n if \"eff_med\" in metrics: \n test_results[\"eff_med\"] = np.median(intervals[:,1]-intervals[:,0])\n if \"time_fit\" in metrics:\n test_results[\"time_fit\"] = self.time_fit\n toc = time.time()\n self.time_evaluate = toc-tic\n if \"time_evaluate\" in metrics:\n test_results[\"time_evaluate\"] = self.time_evaluate\n return test_results\n \nclass ConformalPredictiveSystem(ConformalPredictor):\n \"\"\"\n A conformal predictive system transforms point predictions \n (regression values) into cumulative distribution functions \n (conformal predictive distributions).\n \"\"\"\n \n def __repr__(self):\n if self.fitted:\n return (f\"ConformalPredictiveSystem(fitted={self.fitted}, \"\n f\"normalized={self.normalized}, \"\n f\"mondrian={self.mondrian})\")\n else:\n return f\"ConformalPredictiveSystem(fitted={self.fitted})\"\n\n def fit(self, residuals, sigmas=None, bins=None, seed=None):\n \"\"\"\n Fit conformal predictive system.\n\n Parameters\n ----------\n residuals : array-like of shape (n_values,)\n actual values - predicted values\n sigmas: array-like of shape (n_values,), default=None\n difficulty estimates\n bins : array-like of shape (n_values,), default=None\n Mondrian categories\n seed : int, default=None\n set random seed\n\n Returns\n -------\n self : object\n Fitted ConformalPredictiveSystem.\n\n Examples\n --------\n Assuming that ``y_cal`` and ``y_hat_cal`` are vectors with true and\n predicted targets for some calibration set, then a standard conformal\n predictive system can be formed from the residuals:\n\n .. code-block:: python\n\n residuals_cal = y_cal - y_hat_cal\n\n from crepes import ConformalPredictiveSystem\n\n cps_std = ConformalPredictiveSystem() \n\n cps_std.fit(residuals_cal) \n\n Assuming that ``sigmas_cal`` is a vector with difficulty estimates,\n then a normalized conformal predictive system can be fitted in the \n following way:\n\n .. code-block:: python\n\n cps_norm = ConformalPredictiveSystem()\n cps_norm.fit(residuals_cal, sigmas=sigmas_cal)\n\n Assuming that ``bins_cals`` is a vector with Mondrian categories (bin\n labels), then a Mondrian conformal predictive system can be fitted in\n the following way:\n\n .. code-block:: python\n\n cps_mond = ConformalPredictiveSystem()\n cps_mond.fit(residuals_cal, bins=bins_cal)\n\n A normalized Mondrian conformal predictive system can be fitted in the\n following way:\n\n .. code-block:: python\n\n cps_norm_mond = ConformalPredictiveSystem()\n cps_norm_mond.fit(residuals_cal, sigmas=sigmas_cal, \n bins=bins_cal)\n\n Note\n ----\n By providing a random seed, e.g., ``seed=123``, calls to the methods\n ``predict`` and ``evaluate`` of the :class:`.ConformalPredictiveSystem`\n object will be deterministic. \n \"\"\"\n tic = time.time()\n if bins is None:\n self.mondrian = False\n if sigmas is None:\n self.normalized = False\n self.alphas = np.sort(residuals)\n else:\n self.normalized = True\n self.alphas = np.sort(residuals/sigmas)\n else: \n self.mondrian = True\n bin_values = np.unique(bins)\n if sigmas is None: \n self.normalized = False\n self.alphas = (bin_values, [np.sort(\n residuals[bins==b]) for b in bin_values])\n else:\n self.normalized = True\n self.alphas = (bin_values, [np.sort(\n residuals[bins==b]/sigmas[bins==b]) for b in bin_values]) \n self.fitted = True\n self.seed = seed\n toc = time.time()\n self.time_fit = toc-tic\n return self\n \n def predict(self, y_hat, sigmas=None, bins=None,\n y=None, lower_percentiles=None, higher_percentiles=None,\n y_min=-np.inf, y_max=np.inf, return_cpds=False,\n cpds_by_bins=False, seed=None): \n \"\"\"\n Predict using conformal predictive system.\n\n Parameters\n ----------\n y_hat : array-like of shape (n_values,)\n predicted values\n sigmas : array-like of shape (n_values,), default=None\n difficulty estimates\n bins : array-like of shape (n_values,), default=None\n Mondrian categories\n y : float, int or array-like of shape (n_values,), default=None\n values for which p-values should be returned\n lower_percentiles : array-like of shape (l_values,), default=None\n percentiles for which a lower value will be output \n in case a percentile lies between two values\n (similar to `interpolation=\"lower\"` in `numpy.percentile`)\n higher_percentiles : array-like of shape (h_values,), default=None\n percentiles for which a higher value will be output \n in case a percentile lies between two values\n (similar to `interpolation=\"higher\"` in `numpy.percentile`)\n y_min : float or int, default=-numpy.inf\n The minimum value to include in prediction intervals.\n y_max : float or int, default=numpy.inf\n The maximum value to include in prediction intervals.\n return_cpds : Boolean, default=False\n specifies whether conformal predictive distributions (cpds)\n should be output or not\n cpds_by_bins : Boolean, default=False\n specifies whether the output cpds should be grouped by bin or not; \n only applicable when bins is not None and return_cpds = True\n seed : int, default=None\n set random seed\n\n Returns\n -------\n results : ndarray of shape (n_values, n_cols) or (n_values,)\n the shape is (n_values, n_cols) if n_cols > 1 and otherwise\n (n_values,), where n_cols = p_values+l_values+h_values where \n p_values = 1 if y is not None and 0 otherwise, l_values are the\n number of lower percentiles, and h_values are the number of higher\n percentiles. Only returned if n_cols > 0.\n cpds : ndarray of (n_values, c_values), ndarray of (n_values,)\n or list of ndarrays\n conformal predictive distributions. Only returned if \n return_cpds == True. If bins is None, the distributions are\n represented by a single array, where the number of columns\n (c_values) is determined by the number of residuals of the fitted\n conformal predictive system. Otherwise, the distributions\n are represented by a vector of arrays, if cpds_by_bins = False,\n or a list of arrays, with one element for each bin, if \n cpds_by_bins = True.\n\n Examples\n --------\n Assuming that ``y_hat_test`` and ``y_test`` are vectors with predicted\n and true targets, respectively, for a test set and ``cps_std`` a fitted\n standard conformal predictive system, the p-values for the true targets \n can be obtained by:\n\n .. code-block:: python\n\n p_values = cps_std.predict(y_hat_test, y=y_test)\n\n The p-values with respect to some specific value, e.g., 37, can be\n obtained by:\n\n .. code-block:: python\n\n p_values = cps_std.predict(y_hat_test, y=37)\n\n Assuming that ``sigmas_test`` is a vector with difficulty estimates for\n the test set and ``cps_norm`` a fitted normalized conformal predictive \n system, then the 90th and 95th percentiles can be obtained by:\n\n .. code-block:: python\n\n percentiles = cps_norm.predict(y_hat_test, sigmas=sigmas_test,\n higher_percentiles=[90,95])\n\n In the above example, the nearest higher value is returned, if there is\n no value that corresponds exactly to the requested percentile. If we\n instead would like to retrieve the nearest lower value, we should \n write:\n\n .. code-block:: python\n\n percentiles = cps_norm.predict(y_hat_test, sigmas=sigmas_test,\n lower_percentiles=[90,95])\n\n Assuming that ``bins_test`` is a vector with Mondrian categories (bin \n labels) for the test set and ``cps_mond`` a fitted Mondrian conformal \n regressor, then the following returns prediction intervals at the \n 95% confidence level, where the intervals are lower-bounded by 0:\n\n .. code-block:: python\n\n intervals = cps_mond.predict(y_hat_test, bins=bins_test,\n lower_percentiles=2.5,\n higher_percentiles=97.5,\n y_min=0)\n\n If we would like to obtain the conformal distributions, we could write\n the following:\n\n .. code-block:: python\n\n cpds = cps_norm.predict(y_hat_test, sigmas=sigmas_test,\n return_cpds=True)\n\n The output of the above will be an array with a row for each test\n instance and a column for each calibration instance (residual).\n For a Mondrian conformal predictive system, the above will instead\n result in a vector, in which each element is a vector, as the number\n of calibration instances may vary between categories. If we instead\n would like an array for each category, this can be obtained by:\n\n .. code-block:: python\n\n cpds = cps_norm.predict(y_hat_test, sigmas=sigmas_test,\n return_cpds=True, cpds_by_bins=True)\n\n Note\n ----\n In case the calibration set is too small for the specified lower and\n higher percentiles, a warning will be issued and the output will be \n ``y_min`` and ``y_max``, respectively.\n\n Note\n ----\n Setting ``return_cpds=True`` may consume a lot of memory, as a matrix is\n generated for which the number of elements is the product of the number \n of calibration and test objects, unless a Mondrian approach is employed; \n for the latter, this number is reduced by increasing the number of bins.\n\n Note\n ----\n Setting ``cpds_by_bins=True`` has an effect only for Mondrian conformal \n predictive systems.\n\n Note\n ----\n If a value for ``seed`` is given, it will take precedence over any ``seed``\n value given when calling ``fit``. \n \"\"\"\n tic = time.time()\n if seed is None:\n seed = self.seed\n if seed is not None:\n random_state = np.random.get_state()\n np.random.seed(seed)\n if self.mondrian:\n bin_values, bin_alphas = self.alphas\n bin_indexes = [np.argwhere(bins == b).T[0] for b in bin_values]\n no_prec_result_cols = 0\n if isinstance(lower_percentiles, (int, float, np.integer, np.floating)):\n lower_percentiles = [lower_percentiles]\n if isinstance(higher_percentiles, (int, float, np.integer, np.floating)):\n higher_percentiles = [higher_percentiles]\n if lower_percentiles is None:\n lower_percentiles = []\n if higher_percentiles is None:\n higher_percentiles = []\n if (np.array(lower_percentiles) > 100).any() or \\\n (np.array(lower_percentiles) < 0).any() or \\\n (np.array(higher_percentiles) > 100).any() or \\\n (np.array(higher_percentiles) < 0).any():\n raise ValueError(\"All percentiles must be in the range [0,100]\")\n no_result_columns = \\\n (y is not None) + len(lower_percentiles) + len(higher_percentiles)\n if no_result_columns > 0:\n result = np.zeros((len(y_hat),no_result_columns))\n if y is not None:\n if isinstance(y, pd.Series):\n y = y.values\n no_prec_result_cols += 1\n gammas = np.random.rand(len(y_hat))\n if isinstance(y, (int, float, np.integer, np.floating)):\n if not self.mondrian:\n if self.normalized:\n result[:,0] = np.array(\n [(len(np.argwhere(\n (y_hat[i]+sigmas[i]*self.alphas) 0:\n if not self.mondrian:\n lower_indexes = np.array([int(lower_percentile/100 \\\n * (len(self.alphas)+1))-1\n for lower_percentile in lower_percentiles])\n too_low_indexes = np.argwhere(lower_indexes < 0)\n if len(too_low_indexes) > 0:\n lower_indexes[too_low_indexes[:,0]] = 0\n percentiles_to_show = \" \".join([\n str(lower_percentiles[i])\n for i in too_low_indexes[:,0]])\n warnings.warn(\"the no. of calibration examples is \" \\\n \"too small for the following lower \" \\\n f\"percentiles: {percentiles_to_show}; \"\\\n \"the corresponding values are \" \\\n \"set to y_min\")\n y_min_columns = [no_prec_result_cols+i\n for i in too_low_indexes[:,0]]\n percentile_indexes = lower_indexes\n else:\n too_small_bins = []\n binned_lower_indexes = []\n for b in range(len(bin_values)):\n lower_indexes = np.array([int(lower_percentile/100 \\\n * (len(bin_alphas[b])+1))-1\n for lower_percentile\n in lower_percentiles])\n binned_lower_indexes.append(lower_indexes)\n too_low_indexes = np.argwhere(lower_indexes < 0)\n if len(too_low_indexes) > 0:\n lower_indexes[too_low_indexes[:,0]] = 0\n too_small_bins.append(str(bin_values[b])) \n y_min_columns.append([no_prec_result_cols+i\n for i in too_low_indexes[:,0]])\n else:\n y_min_columns.append([])\n percentile_indexes = [binned_lower_indexes]\n if len(too_small_bins) > 0:\n if len(too_small_bins) < 11:\n bins_to_show = \" \".join(too_small_bins)\n else:\n bins_to_show = \" \".join(\n too_small_bins[:10]+['...'])\n warnings.warn(\"the no. of calibration examples is \" \\\n \"too small for some lower percentile \" \\\n \"in the following bins:\" \\\n f\"{bins_to_show}; \"\\\n \"the corresponding values are \" \\\n \"set to y_min\") \n if len(higher_percentiles) > 0:\n if not self.mondrian:\n higher_indexes = np.array(\n [int(np.ceil(higher_percentile/100 \\\n * (len(self.alphas)+1)))-1\n for higher_percentile in higher_percentiles],\n dtype=int)\n too_high_indexes = np.array(\n [i for i in range(len(higher_indexes))\n if higher_indexes[i] > len(self.alphas)-1], dtype=int)\n if len(too_high_indexes) > 0:\n higher_indexes[too_high_indexes] = len(self.alphas)-1\n percentiles_to_show = \" \".join(\n [str(higher_percentiles[i])\n for i in too_high_indexes])\n warnings.warn(\"the no. of calibration examples is \" \\\n \"too small for the following higher \" \\\n f\"percentiles: {percentiles_to_show}; \"\\\n \"the corresponding values are \" \\\n \"set to y_max\")\n y_max_columns = [no_prec_result_cols+len(lower_percentiles)+i\n for i in too_high_indexes]\n if len(percentile_indexes) == 0:\n percentile_indexes = higher_indexes\n else:\n percentile_indexes = np.concatenate((lower_indexes,\n higher_indexes))\n else:\n too_small_bins = []\n binned_higher_indexes = []\n for b in range(len(bin_values)):\n higher_indexes = np.array([\n int(np.ceil(higher_percentile/100 \\\n * (len(bin_alphas[b])+1)))-1\n for higher_percentile in higher_percentiles])\n binned_higher_indexes.append(higher_indexes)\n too_high_indexes = np.array(\n [i for i in range(len(higher_indexes))\n if higher_indexes[i] > len(bin_alphas[b])-1],\n dtype=int)\n if len(too_high_indexes) > 0:\n higher_indexes[too_high_indexes] = -1\n too_small_bins.append(str(bin_values[b]))\n y_max_columns.append([no_prec_result_cols + \\\n len(lower_percentiles)+i\n for i in too_high_indexes])\n else:\n y_max_columns.append([])\n if len(percentile_indexes) == 0:\n percentile_indexes = [binned_higher_indexes]\n else:\n percentile_indexes.append(binned_higher_indexes)\n if len(too_small_bins) > 0:\n if len(too_small_bins) < 11:\n bins_to_show = \" \".join(too_small_bins)\n else:\n bins_to_show = \" \".join(\n too_small_bins[:10]+['...'])\n warnings.warn(\"the no. of calibration examples is \" \\\n \"too small for some higher percentile \" \\\n \"in the following bins:\" \\\n f\"{bins_to_show}; \"\\\n \"the corresponding values are \" \\\n \"set to y_max\")\n if len(percentile_indexes) > 0:\n if not self.mondrian:\n if self.normalized:\n result[:,no_prec_result_cols:no_prec_result_cols \\\n + len(percentile_indexes)] = np.array(\n [(y_hat[i] + sigmas[i] * \\\n self.alphas)[percentile_indexes]\n for i in range(len(y_hat))])\n else:\n result[:,no_prec_result_cols:no_prec_result_cols \\\n + len(percentile_indexes)] = np.array(\n [(y_hat[i]+self.alphas)[percentile_indexes]\n for i in range(len(y_hat))])\n if len(y_min_columns) > 0:\n result[:,y_min_columns] = y_min\n if len(y_max_columns) > 0:\n result[:,y_max_columns] = y_max\n else:\n if len(percentile_indexes) == 1:\n percentile_indexes = percentile_indexes[0]\n else:\n percentile_indexes = [np.concatenate(\n (percentile_indexes[0][b],percentile_indexes[1][b]))\n for b in range(len(bin_values))]\n if self.normalized:\n for b in range(len(bin_values)):\n if len(bin_indexes[b]) > 0:\n result[bin_indexes[b],\n no_prec_result_cols:no_prec_result_cols \\\n + len(percentile_indexes[b])] = \\\n np.array([(y_hat[i] + sigmas[i] * \\\n bin_alphas[b])[\n percentile_indexes[b]]\n for i in bin_indexes[b]])\n else:\n for b in range(len(bin_values)):\n if len(bin_indexes[b]) > 0:\n result[bin_indexes[b],\n no_prec_result_cols:no_prec_result_cols \\\n + len(percentile_indexes[b])] = np.array(\n [(y_hat[i]+bin_alphas[b])[\n percentile_indexes[b]]\n for i in bin_indexes[b]])\n if len(y_min_columns) > 0:\n for b in range(len(bin_values)):\n if len(bin_indexes[b]) > 0 and \\\n len(y_min_columns[b]) > 0:\n result[bin_indexes[b],y_min_columns[b]] = y_min\n if len(y_max_columns) > 0:\n for b in range(len(bin_values)):\n if len(bin_indexes[b]) > 0 and \\\n len(y_max_columns[b]) > 0:\n result[bin_indexes[b],y_max_columns[b]] = y_max\n if y_min > -np.inf:\n result[:,\n no_prec_result_cols:no_prec_result_cols \\\n + len(percentile_indexes)]\\\n [result[:,no_prec_result_cols:no_prec_result_cols \\\n + len(percentile_indexes)]y_max] = y_max\n no_prec_result_cols += len(percentile_indexes)\n toc = time.time()\n self.time_predict = toc-tic \n if no_result_columns > 0 and result.shape[1] == 1:\n result = result[:,0]\n if return_cpds:\n if not self.mondrian:\n if self.normalized:\n cpds = np.array([y_hat[i]+sigmas[i]*self.alphas\n for i in range(len(y_hat))])\n else:\n cpds = np.array([y_hat[i]+self.alphas\n for i in range(len(y_hat))])\n else: \n if self.normalized:\n cpds = [np.array([y_hat[i]+sigmas[i]*bin_alphas[b]\n for i in bin_indexes[b]])\n for b in range(len(bin_values))]\n else:\n cpds = [np.array([y_hat[i]+bin_alphas[b] for\n i in bin_indexes[b]])\n for b in range(len(bin_values))]\n if no_result_columns > 0 and return_cpds:\n if not self.mondrian or cpds_by_bins:\n cpds_out = cpds\n else:\n cpds_out = np.empty(len(y_hat), dtype=object)\n for b in range(len(bin_values)):\n cpds_out[bin_indexes[b]] = [cpds[b][i]\n for i in range(len(cpds[b]))]\n return result, cpds_out\n elif no_result_columns > 0:\n return result\n elif return_cpds:\n if not self.mondrian or cpds_by_bins:\n cpds_out = cpds\n else:\n cpds_out = np.empty(len(y_hat), dtype=object)\n for b in range(len(bin_values)):\n cpds_out[bin_indexes[b]] = [cpds[b][i]\n for i in range(len(cpds[b]))]\n return cpds_out\n if seed is not None:\n np.random.set_state(random_state)\n\n def evaluate(self, y_hat, y, sigmas=None, bins=None,\n confidence=0.95, y_min=-np.inf, y_max=np.inf,\n metrics=None, seed=None):\n \"\"\"\n Evaluate conformal predictive system.\n\n Parameters\n ----------\n y_hat : array-like of shape (n_values,)\n predicted values\n y : array-like of shape (n_values,)\n correct target values\n sigmas : array-like of shape (n_values,), default=None,\n difficulty estimates\n bins : array-like of shape (n_values,), default=None,\n Mondrian categories\n confidence : float in range (0,1), default=0.95\n confidence level\n y_min : float or int, default=-numpy.inf\n minimum value to include in prediction intervals\n y_max : float or int, default=numpy.inf\n maximum value to include in prediction intervals\n metrics : a string or a list of strings, default=list of all \n metrics; [\"error\", \"eff_mean\",\"eff_med\", \"CRPS\", \"time_fit\",\n \"time_evaluate\"]\n seed : int, default=None\n set random seed\n \n Returns\n -------\n results : dictionary with a key for each selected metric \n estimated performance using the metrics\n\n Examples\n --------\n Assuming that ``y_hat_test`` and ``y_test`` are vectors with predicted\n and true targets for a test set, ``sigmas_test`` and ``bins_test`` are\n vectors with difficulty estimates and Mondrian categories (bin labels) \n for the test set, and ``cps_norm_mond`` is a fitted normalized Mondrian\n conformal predictive system, then the latter can be evaluated at the \n default confidence level with respect to error, mean and median \n efficiency (interval size, given the default confidence level) and \n continuous-ranked probability score (CRPS) by:\n\n .. code-block:: python\n\n results = cps_norm_mond.evaluate(y_hat_test, y_test, \n sigmas=sigmas_test, bins=bins_test,\n metrics=[\"error\", \"eff_mean\", \n \"eff_med\", \"CRPS\"])\n\n Note\n ----\n The use of the metric ``CRPS`` may consume a lot of memory, as a matrix\n is generated for which the number of elements is the product of the \n number of calibration and test objects, unless a Mondrian approach is \n employed; for the latter, this number is reduced by increasing the number \n of bins.\n\n Note\n ----\n If a value for ``seed`` is given, it will take precedence over any ``seed``\n value given when calling ``fit``. \n \"\"\"\n tic = time.time()\n if seed is None:\n seed = self.seed\n if seed is not None:\n random_state = np.random.get_state()\n np.random.seed(seed)\n if isinstance(y, pd.Series):\n y = y.values\n if isinstance(y_hat, pd.Series):\n y_hat = y_hat.values\n test_results = {}\n lower_percentile = (1-confidence)/2*100\n higher_percentile = (confidence+(1-confidence)/2)*100\n if metrics is None:\n metrics = [\"error\",\"eff_mean\",\"eff_med\",\"CRPS\",\"time_fit\",\n \"time_evaluate\"]\n if \"CRPS\" in metrics:\n results, cpds = self.predict(y_hat, sigmas=sigmas, bins=bins, y=y,\n lower_percentiles=lower_percentile,\n higher_percentiles=higher_percentile,\n y_min=y_min, y_max=y_max,\n return_cpds=True, cpds_by_bins=True)\n intervals = results[:,[1,2]]\n if not self.mondrian:\n if self.normalized:\n crps = calculate_crps(cpds, self.alphas, sigmas, y)\n else:\n crps = calculate_crps(cpds, self.alphas,\n np.ones(len(y_hat)), y)\n else:\n bin_values, bin_alphas = self.alphas\n bin_indexes = [np.argwhere(bins == b).T[0]\n for b in bin_values]\n if self.normalized:\n crps = np.sum([calculate_crps(cpds[b],\n bin_alphas[b],\n sigmas[bin_indexes[b]],\n y[bin_indexes[b]]) \\\n * len(bin_indexes[b])\n for b in range(len(bin_values))])/len(y)\n else:\n crps = np.sum([calculate_crps(cpds[b],\n bin_alphas[b],\n np.ones(len(bin_indexes[b])),\n y[bin_indexes[b]]) \\\n * len(bin_indexes[b])\n for b in range(len(bin_values))])/len(y)\n else:\n intervals = self.predict(y_hat, sigmas=sigmas, bins=bins,\n lower_percentiles=lower_percentile,\n higher_percentiles=higher_percentile,\n y_min=y_min, y_max=y_max,\n return_cpds=False)\n if \"error\" in metrics:\n test_results[\"error\"] = 1-np.mean(np.logical_and(\n intervals[:,0]<=y, y<=intervals[:,1]))\n if \"eff_mean\" in metrics: \n test_results[\"eff_mean\"] = np.mean(intervals[:,1]-intervals[:,0])\n if \"eff_med\" in metrics: \n test_results[\"eff_med\"] = np.median(intervals[:,1]-intervals[:,0])\n if \"CRPS\" in metrics:\n test_results[\"CRPS\"] = crps\n if \"time_fit\" in metrics:\n test_results[\"time_fit\"] = self.time_fit\n toc = time.time()\n if seed is not None:\n np.random.set_state(random_state)\n self.time_evaluate = toc-tic\n if \"time_evaluate\" in metrics:\n test_results[\"time_evaluate\"] = self.time_evaluate\n return test_results\n \ndef calculate_crps(cpds, alphas, sigmas, y):\n \"\"\"\n Calculate mean continuous-ranked probability score (crps)\n for a set of conformal predictive distributions.\n\n Parameters\n ----------\n cpds : array-like of shape (n_values, c_values)\n conformal predictive distributions\n alphas : array-like of shape (c_values,)\n sorted (normalized) residuals of the calibration examples \n sigmas : array-like of shape (n_values,),\n difficulty estimates\n y : array-like of shape (n_values,)\n correct target values\n \n Returns\n -------\n crps : float\n mean continuous-ranked probability score for the conformal\n predictive distributions \n \"\"\"\n if len(cpds) > 0:\n widths = np.array([alphas[i+1]-alphas[i] for i in range(len(alphas)-1)])\n cum_probs = np.cumsum([1/len(alphas) for i in range(len(alphas)-1)])\n lower_errors = cum_probs**2\n higher_errors = (1-cum_probs)**2\n cpd_indexes = [np.argwhere(cpds[i] 1 and otherwise\n (n_samples,), where n_cols = p_values+l_values+h_values where \n p_values = 1 if y is not None and 0 otherwise, l_values are the\n number of lower percentiles, and h_values are the number of higher\n percentiles. Only returned if n_cols > 0.\n cpds : ndarray of (n_samples, c_values), ndarray of (n_samples,)\n or list of ndarrays\n conformal predictive distributions. Only returned if \n return_cpds == True. For non-Mondrian conformal predictive systems,\n the distributions are represented by a single array, where the \n number of columns (c_values) is determined by the number of \n residuals of the fitted conformal predictive system. For Mondrian\n conformal predictive systems, the distributions are represented by\n a vector of arrays, if cpds_by_bins = False, or a list of arrays, \n with one element for each Mondrian category, if cpds_by_bins = True.\n\n Examples\n --------\n Assuming that ``X_test`` is a set of test objects, ``y_test`` is a \n vector with true targets, ``w`` is a :class:`.WrapRegressor` object \n calibrated with the option ``cps=True``, p-values for the true targets \n can be obtained by:\n\n .. code-block:: python\n\n p_values = w.predict_cps(X_test, y=y_test)\n\n P-values with respect to some specific value, e.g., 37, can be\n obtained by:\n\n .. code-block:: python\n\n p_values = w.predict_cps(X_test, y=37)\n\n The 90th and 95th percentiles can be obtained by:\n\n .. code-block:: python\n\n percentiles = w.predict_cps(X_test, higher_percentiles=[90,95])\n\n In the above example, the nearest higher value is returned, if there is\n no value that corresponds exactly to the requested percentile. If we\n instead would like to retrieve the nearest lower value, we should \n write:\n\n .. code-block:: python\n\n percentiles = w.predict_cps(X_test, lower_percentiles=[90,95])\n\n The following returns prediction intervals at the 95% confidence level,\n where the intervals are lower-bounded by 0:\n\n .. code-block:: python\n\n intervals = w.predict_cps(X_test,\n lower_percentiles=2.5,\n higher_percentiles=97.5,\n y_min=0)\n\n If we would like to obtain the conformal distributions, we could write\n the following:\n\n .. code-block:: python\n\n cpds = w.predict_cps(X_test, return_cpds=True)\n\n The output of the above will be an array with a row for each test\n instance and a column for each calibration instance (residual).\n If the learner is wrapped with a Mondrian conformal predictive system, \n the above will instead result in a vector, in which each element is a\n vector, as the number of calibration instances may vary between \n categories. If we instead would like an array for each category, this \n can be obtained by:\n\n .. code-block:: python\n\n cpds = w.predict_cps(X_test, return_cpds=True, cpds_by_bins=True)\n\n Note\n ----\n This method is available only if the learner has been wrapped with a\n :class:`.ConformalPredictiveSystem`, i.e., :meth:`.calibrate`\n has been called with the option ``cps=True``.\n\n Note\n ----\n In case the calibration set is too small for the specified lower and\n higher percentiles, a warning will be issued and the output will be \n ``y_min`` and ``y_max``, respectively.\n\n Note\n ----\n Setting ``return_cpds=True`` may consume a lot of memory, as a matrix is\n generated for which the number of elements is the product of the number \n of calibration and test objects, unless a Mondrian approach is employed; \n for the latter, this number is reduced by increasing the number of bins.\n\n Note\n ----\n Setting ``cpds_by_bins=True`` has an effect only for Mondrian conformal \n predictive systems.\n\n Note\n ----\n If a value for ``seed`` is given, it will take precedence over any ``seed``\n value given when calling ``calibrate``. \n \"\"\"\n if isinstance(y, pd.Series):\n y = y.values\n if self.cps is None:\n raise RuntimeError((\"predict_cps requires that calibrate has been\"\n \"called first with cps=True\"))\n else:\n if seed is None:\n seed = self.seed\n if seed is not None:\n random_state = np.random.get_state()\n np.random.seed(seed)\n y_hat = self.learner.predict(X)\n if self.de is None:\n sigmas = None\n else:\n sigmas = self.de.apply(X)\n if self.mc is None:\n bins = None\n elif isinstance(self.mc, MondrianCategorizer):\n bins = self.mc.apply(X)\n else:\n bins = self.mc(X)\n result = self.cps.predict(y_hat, sigmas=sigmas, bins=bins,\n y=y, lower_percentiles=lower_percentiles,\n higher_percentiles=higher_percentiles,\n y_min=y_min, y_max=y_max,\n return_cpds=return_cpds,\n cpds_by_bins=cpds_by_bins)\n if seed is not None:\n np.random.set_state(random_state)\n return result\n\n def evaluate(self, X, y, confidence=0.95, y_min=-np.inf, y_max=np.inf,\n metrics=None, seed=None):\n \"\"\"\n Evaluate :class:`.ConformalRegressor` or \n :class:`.ConformalPredictiveSystem`.\n\n Parameters\n ----------\n X : array-like of shape (n_samples, n_features)\n set of objects\n y : array-like of shape (n_samples,)\n correct target values\n confidence : float in range (0,1), default=0.95\n confidence level\n y_min : float or int, default=-numpy.inf\n minimum value to include in prediction intervals\n y_max : float or int, default=numpy.inf\n maximum value to include in prediction intervals\n metrics : a string or a list of strings, default=list of all \n metrics; for a learner wrapped with a conformal regressor\n these are \"error\", \"eff_mean\",\"eff_med\", \"time_fit\", and\n \"time_evaluate\", while if wrapped with a conformal predictive\n system, the metrics also include \"CRPS\". \n seed : int, default=None\n set random seed\n \n Returns\n -------\n results : dictionary with a key for each selected metric \n estimated performance using the metrics \n\n Examples\n --------\n Assuming that ``X_test`` is a set of test objects, ``y_test`` is a \n vector with true targets, and ``w`` is a calibrated \n :class:`.WrapRegressor` object, then the latter can be evaluated at \n the 90% confidence level with respect to error, mean and median \n efficiency (interval size) by:\n\n .. code-block:: python\n\n results = w.evaluate(X_test, y_test, confidence=0.9,\n metrics=[\"error\", \"eff_mean\", \"eff_med\"])\n\n Note\n ----\n If included in the list of metrics, \"CRPS\" (continuous-ranked\n probability score) will be ignored if the :class:`.WrapRegressor` object\n has been calibrated with the (default) option ``cps=False``, i.e., the \n learner is wrapped with a :class:`.ConformalRegressor`.\n\n Note\n ----\n The use of the metric ``CRPS`` may consume a lot of memory, as a matrix\n is generated for which the number of elements is the product of the \n number of calibration and test objects, unless a Mondrian approach is \n employed; for the latter, this number is reduced by increasing the number \n of categories.\n\n Note\n ----\n The reported result for ``time_fit`` only considers fitting the\n conformal regressor or predictive system; not for fitting the\n learner.\n\n Note\n ----\n If a value for ``seed`` is given, it will take precedence over any ``seed``\n value given when calling ``calibrate``.\n \"\"\"\n if isinstance(y, pd.Series):\n y = y.values\n if not self.calibrated:\n raise RuntimeError((\"evaluate requires that calibrate has been\"\n \"called first\"))\n else:\n if seed is None:\n seed = self.seed\n if seed is not None:\n random_state = np.random.get_state()\n np.random.seed(seed)\n y_hat = self.learner.predict(X)\n if self.de is None:\n sigmas = None\n else:\n sigmas = self.de.apply(X)\n if self.mc is None:\n bins = None\n elif isinstance(self.mc, MondrianCategorizer):\n bins = self.mc.apply(X)\n else:\n bins = self.mc(X)\n if self.cr is not None:\n result = self.cr.evaluate(y_hat, y, sigmas=sigmas,\n bins=bins, confidence=confidence,\n y_min=y_min, y_max=y_max)\n else:\n result = self.cps.evaluate(y_hat, y, sigmas=sigmas,\n bins=bins, confidence=confidence,\n y_min=y_min, y_max=y_max)\n if seed is not None:\n np.random.set_state(random_state)\n return result","metadata":{"trusted":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2024-11-06T01:51:26.677886Z","iopub.execute_input":"2024-11-06T01:51:26.678467Z","iopub.status.idle":"2024-11-06T01:51:27.012953Z","shell.execute_reply.started":"2024-11-06T01:51:26.678405Z","shell.execute_reply":"2024-11-06T01:51:27.011680Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class MyLogger:\n \"\"\"\n This class helps to suppress logs in lightgbm and Optuna\n Source - https://github.com/microsoft/LightGBM/issues/6014\n \"\"\"\n\n def init(self, logging_lbl: str):\n self.logger = logging.getLogger(logging_lbl)\n self.logger.setLevel(logging.ERROR)\n\n def info(self, message):\n pass\n\n def warning(self, message):\n pass\n\n def error(self, message):\n self.logger.error(message)\n\nl = MyLogger()\nl.init(logging_lbl = \"lightgbm_custom\")\nlgb.register_logger(l)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.015263Z","iopub.execute_input":"2024-11-06T01:51:27.015627Z","iopub.status.idle":"2024-11-06T01:51:27.024412Z","shell.execute_reply.started":"2024-11-06T01:51:27.015589Z","shell.execute_reply":"2024-11-06T01:51:27.022907Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"season_dtype = pl.Enum(['Spring', 'Summer', 'Fall', 'Winter','Missing'])\ntarget_labels = ['None', 'Mild', 'Moderate', 'Severe']\n\ntrain = (\n pl.read_csv('/kaggle/input/child-mind-institute-problematic-internet-use/train.csv')\n .with_columns(\n pl.col('^.*Season$').fill_null('Missing').cast(season_dtype),\n )\n # .drop('^PCIAT.*$','^PAQ_*$')\n)\n\ntest = (\n pl.read_csv('/kaggle/input/child-mind-institute-problematic-internet-use/test.csv')\n .with_columns(pl.col('^.*Season$').fill_null('Missing').cast(season_dtype))\n .drop('^PCIAT.*$','^PAQ_*$')\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.026046Z","iopub.execute_input":"2024-11-06T01:51:27.026541Z","iopub.status.idle":"2024-11-06T01:51:27.054631Z","shell.execute_reply.started":"2024-11-06T01:51:27.026483Z","shell.execute_reply":"2024-11-06T01:51:27.053423Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def recalculate_sii(row):\n if pd.isna(row['PCIAT-PCIAT_Total']):\n return np.nan\n max_possible = row['PCIAT-PCIAT_Total'] + row[PCIAT_cols].isna().sum() * 5\n if row['PCIAT-PCIAT_Total'] <= 30 and max_possible <= 30:\n return 0\n elif 31 <= row['PCIAT-PCIAT_Total'] <= 49 and max_possible <= 49:\n return 1\n elif 50 <= row['PCIAT-PCIAT_Total'] <= 79 and max_possible <= 79:\n return 2\n elif row['PCIAT-PCIAT_Total'] >= 80 and max_possible >= 80:\n return 3\n return np.nan\n\nif FIX_SII:\n PCIAT_cols = [f'PCIAT-PCIAT_{i+1:02d}' for i in range(20)]\n \n train = train.to_pandas()\n train['sii'] = (\n train.apply(recalculate_sii, axis=1)\n \n )\n\nif isinstance(train, pd.DataFrame):\n train = pl.from_pandas(train).drop('^PCIAT.*$','^PAQ_*$')\nelse:\n train = train.drop('^PCIAT.*$','^PAQ_*$')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.056464Z","iopub.execute_input":"2024-11-06T01:51:27.057780Z","iopub.status.idle":"2024-11-06T01:51:27.070625Z","shell.execute_reply.started":"2024-11-06T01:51:27.057713Z","shell.execute_reply":"2024-11-06T01:51:27.069287Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"missing = train.select([(pl.col(c).is_null().mean()).alias(c) for c in train.columns])\n\nmissing = (\n missing.transpose(include_header=True) \n .filter(pl.col(\"column_0\") < 0.5) \n .with_columns((pl.col(\"column_0\") * 100).round(4).alias('%')) \n .sort('%', descending=True) \n)\n\nprint(f\"{missing.height} columns have less than 50% missing.\")","metadata":{"execution":{"iopub.status.busy":"2024-11-06T01:51:27.074502Z","iopub.execute_input":"2024-11-06T01:51:27.074994Z","iopub.status.idle":"2024-11-06T01:51:27.088667Z","shell.execute_reply.started":"2024-11-06T01:51:27.074942Z","shell.execute_reply":"2024-11-06T01:51:27.087063Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cols = [i for i in train.columns if i in missing['column'].to_numpy()]\ndftr = train[cols]\ndfte = test[cols[:-1]]\n\n# DON'T WANT TO USE SII.. yet anyways.\ndftr = dftr.filter(~pl.col('sii').is_null())","metadata":{"execution":{"iopub.status.busy":"2024-11-06T01:51:27.092731Z","iopub.execute_input":"2024-11-06T01:51:27.093189Z","iopub.status.idle":"2024-11-06T01:51:27.102684Z","shell.execute_reply.started":"2024-11-06T01:51:27.093145Z","shell.execute_reply":"2024-11-06T01:51:27.101257Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def conformal_intervals(X, y, X_test, n_splits=5, confidence=0.90, seed=SEED):\n kf = KFold(n_splits=n_splits, shuffle=True, random_state=seed)\n intervals = []\n \n for train_index, cal_index in kf.split(X):\n X_prop_train, X_cal = X[train_index], X[cal_index]\n y_prop_train, y_cal = y[train_index], y[cal_index]\n \n train_pool = Pool(X_prop_train.astype(np.float32), y_prop_train)\n cal_pool = Pool(X_cal.astype(np.float32), y_cal)\n \n model = CatBoostRegressor(\n iterations=10000,\n loss_function='RMSE',\n depth=4,\n subsample=0.8,\n colsample_bylevel=0.8,\n random_seed=seed,\n verbose=0\n )\n \n model.fit(train_pool, eval_set=cal_pool, early_stopping_rounds=30)\n \n m = WrapRegressor(model)\n m.calibrate(X_cal, y_cal)\n \n interval = m.predict_int(X_test, confidence=confidence)\n intervals.append(interval)\n \n return np.mean(intervals, axis=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.104842Z","iopub.execute_input":"2024-11-06T01:51:27.105610Z","iopub.status.idle":"2024-11-06T01:51:27.117567Z","shell.execute_reply.started":"2024-11-06T01:51:27.105545Z","shell.execute_reply":"2024-11-06T01:51:27.116294Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def feature_eng(dftr, features):\n median_df = (\n dftr.group_by(\"Basic_Demos-Age\")\n .agg(\n (pl.col(\"Physical-Height\").median()).alias(\"Feat_Height\"),\n (pl.col(\"Physical-Weight\").median()).alias(\"Feat_Weight\"),\n )\n )\n\n dftr = dftr.join(median_df, on=\"Basic_Demos-Age\", how=\"left\").with_columns(\n pl.when(pl.col(\"Physical-Height\").is_null())\n .then(pl.col(\"Feat_Height\"))\n .otherwise(pl.col(\"Physical-Height\"))\n .alias(\"Physical-Height\"),\n \n pl.when(pl.col(\"Physical-Weight\").is_null())\n .then(pl.col(\"Feat_Weight\"))\n .otherwise(pl.col(\"Physical-Weight\"))\n .alias(\"Physical-Weight\"),\n ).drop([\"Feat_Height\", \"Feat_Weight\"])\n\n dftr = dftr.with_columns([\n pl.when(pl.col('Basic_Demos-Age') < 10)\n .then(0)\n .when((pl.col('Basic_Demos-Age') >= 10) & (pl.col('Basic_Demos-Age') < 14))\n .then(1)\n .otherwise(2)\n .alias('age_group')\n ])\n\n # for FEAT in features:\n # tmp = dftr.filter(~pl.col(FEAT).is_null())\n # X = tmp.drop(FEAT).to_numpy()\n # y = tmp[FEAT].to_numpy()\n \n # X_test = dftr.drop(FEAT).to_numpy()\n # pred_int = conformal_intervals(X, y, X_test)\n \n # dftr = dftr.with_columns([\n # pl.lit(pred_int[:, 0]).alias(f'{FEAT}_0'),\n # pl.lit(pred_int[:, 1]).alias(f'{FEAT}_1')\n # ])\n \n return dftr","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.135943Z","iopub.execute_input":"2024-11-06T01:51:27.136507Z","iopub.status.idle":"2024-11-06T01:51:27.149965Z","shell.execute_reply.started":"2024-11-06T01:51:27.136420Z","shell.execute_reply":"2024-11-06T01:51:27.148625Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"unique_ids = dftr.to_pandas().drop_duplicates('id')[['id', 'sii']]\nIDS = sorted(list(unique_ids['id']))\n\nFOLDS = []\nstrat_kf = StratifiedKFold(\n n_splits=N_FOLDS,\n shuffle=True,\n random_state=SEED\n)\n\ninitial_fold_ids = []\nfor _, (train_idx, val_idx) in enumerate(strat_kf.split(unique_ids['id'], unique_ids['sii'])):\n val_ids = unique_ids['id'].iloc[val_idx].tolist()\n initial_fold_ids.append(val_ids)\n\nFOLDS.append(initial_fold_ids)\n\nnp.random.seed(SEED)\nfor _ in range(N_BAGS - 1):\n random_fold_assignment = np.random.randint(0, N_FOLDS, len(IDS))\n bag_folds = [np.array(IDS)[random_fold_assignment == fold].tolist() for fold in range(N_FOLDS)]\n FOLDS.append(bag_folds)\n\n# ------------ prepare\ntarget='sii'\n\ny = dftr.get_column(target)\nX = dftr.drop(['id', target]).to_pandas()\n\na, b = y.mean(), y.var(ddof=0)\ny_min, y_max = y.min(), y.max()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.151530Z","iopub.execute_input":"2024-11-06T01:51:27.151986Z","iopub.status.idle":"2024-11-06T01:51:27.234854Z","shell.execute_reply.started":"2024-11-06T01:51:27.151940Z","shell.execute_reply":"2024-11-06T01:51:27.233604Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# dftr.to_pandas().isnull().sum() / len(dftr)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.236709Z","iopub.execute_input":"2024-11-06T01:51:27.237140Z","iopub.status.idle":"2024-11-06T01:51:27.243877Z","shell.execute_reply.started":"2024-11-06T01:51:27.237098Z","shell.execute_reply":"2024-11-06T01:51:27.241169Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def cross_validate_model(\n bst,\n model_train_func,\n n_seeds=N_SEEDS,\n clip_range=(0, 3),\n init_score=2.0,\n):\n global dftr, dfte\n\n gains = np.zeros(len(bst) + 1)\n split = np.zeros(len(bst) + 1)\n feature_names = bst + ['age_group']\n \n tr_preds = [[] for _ in range(len(X))]\n oof_raw_preds = [[] for _ in range(len(y))]\n te_preds = [[] for _ in range(len(dfte))]\n n_splits = N_BAGS * N_FOLDS * n_seeds \n diff = 0\n \n for bag_idx, bag in enumerate(FOLDS):\n for fold_idx in range(N_FOLDS):\n fold_num = bag_idx * N_FOLDS + fold_idx\n\n valid_ids = bag[fold_idx]\n train_ids = [i for fold in bag if fold != valid_ids for i in fold]\n\n idx_tr = dftr.to_pandas().index[dftr.to_pandas()['id'].isin(train_ids)].to_numpy()\n idx_va = dftr.to_pandas().index[dftr.to_pandas()['id'].isin(valid_ids)].to_numpy()\n\n X_tr, X_va = X[bst].iloc[idx_tr], X[bst].iloc[idx_va]\n y_tr, y_va = y[idx_tr], y[idx_va]\n\n # ------------- Feature Engineering\n X_tr, X_va = (\n feature_eng(pl.from_pandas(X_tr), _).to_pandas(),\n feature_eng(pl.from_pandas(X_va), _).to_pandas()\n )\n \n dfte_tmp = feature_eng(dfte[bst], _)\n X_te = dfte_tmp.to_pandas()\n\n for seed in range(n_seeds):\n model = model_train_func(\n X_tr, y_tr, X_va, y_va, init_score, seed=seed\n )\n \n # ------------- Training Predictions\n tr_pred = np.clip(model.predict(X_tr) + init_score, *clip_range)\n tr_pred_rounded = np.round(tr_pred).astype(int)\n for idx, pred in zip(idx_tr, tr_pred_rounded):\n tr_preds[idx].append(pred)\n\n # ------------- Validation Predictions\n y_pred = np.clip(model.predict(X_va) + init_score, *clip_range)\n y_pred_rounded = np.round(y_pred).astype(int)\n for idx, pred in zip(idx_va, y_pred_rounded):\n oof_raw_preds[idx].append(pred)\n\n # ------------- Test Predictions\n te_pred = np.clip(model.predict(X_te) + init_score, *clip_range)\n te_pred_rounded = np.round(te_pred).astype(int)\n for idx, pred in enumerate(te_pred_rounded):\n te_preds[idx].append(pred)\n \n # ------------- Feature Importance\n if hasattr(model, \"feature_importance\"):\n gains += model.feature_importance(importance_type='gain')\n split += model.feature_importance(importance_type='split')\n \n if hasattr(model, \"get_feature_importance\"):\n feature_importances = model.get_feature_importance()\n gains += np.array(feature_importances)\n split += np.array(feature_importances) \n\n trn_eval = cohen_kappa_score(y_tr, tr_pred_rounded, weights='quadratic')\n val_eval = cohen_kappa_score(y_va, y_pred_rounded, weights='quadratic')\n\n diff += abs(trn_eval - val_eval)\n \n total_fold_num = (bag_idx * N_FOLDS * n_seeds) + (fold_idx * n_seeds) + seed\n print(f\"# Processing Fold {total_fold_num + 1}/{n_splits}...\")\n print(f\"TRAIN: {trn_eval:.3f}, EVAL: {val_eval:.3f}, DIFF: {abs(trn_eval - val_eval):.3f}\")\n print()\n\n tr_preds_majority = np.zeros(len(X), dtype=int)\n for idx, preds in enumerate(tr_preds):\n if preds: \n tr_preds_majority[idx] = mode(preds).mode\n\n oof_preds_majority = np.zeros(len(y), dtype=int)\n for idx, preds in enumerate(oof_raw_preds):\n if preds:\n oof_preds_majority[idx] = mode(preds).mode\n\n te_preds_majority = np.zeros(len(dfte), dtype=int)\n for idx, preds in enumerate(te_preds):\n if preds:\n te_preds_majority[idx] = mode(preds).mode\n\n oof_score = cohen_kappa_score(y, oof_preds_majority, weights='quadratic')\n tr_score = cohen_kappa_score(y, tr_preds_majority, weights='quadratic')\n\n print('~' * 50)\n print(f\"# TRAIN: {tr_score:.3f}\")\n print(f\"# OOF: {oof_score:.3f}\")\n print(f\"# DIFF: {(diff / n_splits):.3f}\")\n print('~' * 50)\n\n gains_avg = gains / n_splits \n split_avg = split / n_splits \n\n feature_importance = pd.DataFrame({\n 'feature': feature_names,\n 'gains': gains_avg,\n 'split': split_avg\n })\n\n feature_importance = feature_importance.sort_values(\n by='gains',\n ascending=False\n ).reset_index(drop=True)\n\n return {\n \"oof\": oof_preds_majority,\n \"trn\": tr_preds_majority,\n \"test\": te_preds_majority,\n \"oof_score\": oof_score,\n \"train_score\": tr_score,\n \"feature_importance\": feature_importance,\n }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.247638Z","iopub.execute_input":"2024-11-06T01:51:27.248136Z","iopub.status.idle":"2024-11-06T01:51:27.280859Z","shell.execute_reply.started":"2024-11-06T01:51:27.248094Z","shell.execute_reply":"2024-11-06T01:51:27.279518Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def quadratic_weighted_kappa(preds, data):\n y_true = data.get_label()\n y_pred = preds.clip(y_min, y_max).round()\n qwk = cohen_kappa_score(y_true, y_pred, weights=\"quadratic\")\n return 'QWK', qwk, True\n \ndef qwk_obj(preds, dtrain):\n labels = dtrain.get_label()\n preds = preds.clip(y_min, y_max)\n f = 1/2 * np.sum((preds - labels)**2)\n g = 1/2 * np.sum((preds - a)**2 + b)\n df = preds - labels\n dg = preds - a\n grad = (df/g - f*dg/g**2)*len(labels)\n hess = np.ones(len(labels))\n return grad, hess\n\ndef train_lgb_model(X_tr, y_tr, X_va, y_va, init_score, seed=0):\n lgb_train = lgb.Dataset(X_tr, label=y_tr.to_numpy(),\n init_score=[init_score]*len(X_tr))\n \n lgb_val = lgb.Dataset(X_va, label=y_va.to_numpy(),\n init_score=[init_score]*len(X_va))\n\n p1 = {\n 'lambda_l1': 1.6720583732596097,\n 'lambda_l2': 0.7364155944505759,\n 'feature_fraction': 0.4545944254042284,\n 'colsample_bytree': 0.148464992140102,\n 'max_bin': 55,\n 'max_depth': 4,\n 'num_leaves': 8,\n 'min_data_in_leaf': 77,\n 'path_smooth': 2.4481868043340818,\n 'objective':qwk_obj,\n 'metric': None,\n 'random_state':seed,\n 'extra_trees':True\n }\n\n p2 = dict(\n objective=qwk_obj,\n metric=\"None\",\n feature_fraction=0.5,\n colsample_bytree=0.2,\n max_bin=124,\n max_depth=5,\n num_leaves=2 ** 4,\n min_data_in_leaf=50,\n path_smooth=4,\n verbosity=-1,\n random_state=seed,\n extra_trees=True\n )\n\n model = lgb.train(\n p1,\n lgb_train,\n valid_sets=[lgb_val],\n num_boost_round=10000,\n feval=quadratic_weighted_kappa,\n callbacks=[\n lgb.early_stopping(100, verbose=-1),\n lgb.log_evaluation(0)\n ]\n )\n return model\n \nbst = [\n 'PreInt_EduHx-computerinternet_hoursday',\n 'SDS-SDS_Total_T',\n 'SDS-SDS_Total_Raw',\n 'Basic_Demos-Sex',\n 'Basic_Demos-Age',\n 'Physical-Weight',\n 'Physical-Height',\n 'FGC-FGC_CU',\n 'FGC-FGC_PU',\n 'FGC-FGC_SRL_Zone',\n]\n\nlgb_results = cross_validate_model(bst, train_lgb_model)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:51:27.284106Z","iopub.execute_input":"2024-11-06T01:51:27.285260Z","iopub.status.idle":"2024-11-06T01:52:08.430866Z","shell.execute_reply.started":"2024-11-06T01:51:27.285195Z","shell.execute_reply":"2024-11-06T01:52:08.429701Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# def quadratic_weighted_kappa(preds, data):\n# y_true = data.get_label()\n# y_pred = preds.clip(y_min, y_max).round()\n# qwk = cohen_kappa_score(y_true, y_pred, weights=\"quadratic\")\n# return 'QWK', qwk, True\n \n# def qwk_obj(preds, dtrain):\n# labels = dtrain.get_label()\n# preds = preds.clip(y_min, y_max)\n# f = 1/2 * np.sum((preds - labels)**2)\n# g = 1/2 * np.sum((preds - a)**2 + b)\n# df = preds - labels\n# dg = preds - a\n# grad = (df/g - f*dg/g**2)*len(labels)\n# hess = np.ones(len(labels))\n# return grad, hess\n\n# def tune_train_lgb_model(X_tr, y_tr, X_va, y_va, init_score, p, seed=0):\n# lgb_train = lgb.Dataset(X_tr, label=y_tr.to_numpy(),\n# init_score=[init_score]*len(X_tr))\n \n# lgb_val = lgb.Dataset(X_va, label=y_va.to_numpy(),\n# init_score=[init_score]*len(X_va))\n\n\n\n# model = lgb.train(\n# p,\n# lgb_train,\n# valid_sets=[lgb_val],\n# num_boost_round=10000,\n# feval=quadratic_weighted_kappa,\n# callbacks=[\n# lgb.early_stopping(100, verbose=-1),\n# lgb.log_evaluation(0)\n# ]\n# )\n# return model\n \n# bst = [\n# 'PreInt_EduHx-computerinternet_hoursday',\n# 'SDS-SDS_Total_T',\n# 'SDS-SDS_Total_Raw',\n# 'Basic_Demos-Sex',\n# 'Basic_Demos-Age',\n# 'Physical-Weight',\n# 'Physical-Height',\n# 'FGC-FGC_CU',\n# 'FGC-FGC_PU',\n# 'FGC-FGC_SRL_Zone',\n# ]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.436339Z","iopub.execute_input":"2024-11-06T01:52:08.436707Z","iopub.status.idle":"2024-11-06T01:52:08.443210Z","shell.execute_reply.started":"2024-11-06T01:52:08.436667Z","shell.execute_reply":"2024-11-06T01:52:08.442061Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import optuna\n# from optuna import Trial\n# from scipy.stats import mode\n# import numpy as np\n# import pandas as pd\n\n# def tune_cross_validate_model(trial: Trial, bst, init_score=2.0):\n# params = dict(\n# objective=qwk_obj, \n# metric=\"None\", \n# lambda_l1=trial.suggest_float(\"lambda_l1\", 0.5, 2.0), \n# lambda_l2=trial.suggest_float(\"lambda_l2\", 0.7, 0.8), \n# feature_fraction=trial.suggest_float(\"feature_fraction\", 0.45, 0.55), \n# colsample_bytree=trial.suggest_float(\"colsample_bytree\", 0.1, 0.2), \n# max_bin=trial.suggest_int(\"max_bin\", 50, 250), \n# max_depth=trial.suggest_int(\"max_depth\", 4, 4), \n# num_leaves=trial.suggest_int(\"num_leaves\", 8, 20), \n# min_data_in_leaf=trial.suggest_int(\"min_data_in_leaf\", 50, 80), \n# path_smooth=trial.suggest_float(\"path_smooth\", 2, 4), \n# verbosity=-1,\n# extra_trees=True\n# )\n\n# tr_preds = [[] for _ in range(len(X))]\n# oof_raw_preds = [[] for _ in range(len(y))]\n# te_preds = [[] for _ in range(len(dfte))]\n# n_splits = N_BAGS * N_FOLDS * N_SEEDS \n# diff = 0\n# gains = np.zeros(len(bst) + 1)\n# split = np.zeros(len(bst) + 1)\n# feature_names = bst + ['age_group']\n\n# for bag_idx, bag in enumerate(FOLDS):\n# for fold_idx in range(N_FOLDS):\n# valid_ids = bag[fold_idx]\n# train_ids = [i for fold in bag if fold != valid_ids for i in fold]\n\n# idx_tr = dftr.to_pandas().index[dftr.to_pandas()['id'].isin(train_ids)].to_numpy()\n# idx_va = dftr.to_pandas().index[dftr.to_pandas()['id'].isin(valid_ids)].to_numpy()\n\n# X_tr, X_va = X[bst].iloc[idx_tr], X[bst].iloc[idx_va]\n# y_tr, y_va = y[idx_tr], y[idx_va]\n\n# # ------------- Feature Engineering\n# X_tr, X_va = (\n# feature_eng(pl.from_pandas(X_tr), _).to_pandas(),\n# feature_eng(pl.from_pandas(X_va), _).to_pandas()\n# )\n# dfte_tmp = feature_eng(dfte[bst], _)\n# X_te = dfte_tmp.to_pandas()\n\n# for seed in range(N_SEEDS):\n# model = tune_train_lgb_model(X_tr, y_tr, X_va, y_va, init_score, params, seed=seed)\n \n# # ------------- Training Predictions\n# tr_pred = np.clip(model.predict(X_tr) + init_score, 0, 3)\n# tr_pred_rounded = np.round(tr_pred).astype(int)\n# for idx, pred in zip(idx_tr, tr_pred_rounded):\n# tr_preds[idx].append(pred)\n\n# # ------------- Validation Predictions\n# y_pred = np.clip(model.predict(X_va) + init_score, 0, 3)\n# y_pred_rounded = np.round(y_pred).astype(int)\n# for idx, pred in zip(idx_va, y_pred_rounded):\n# oof_raw_preds[idx].append(pred)\n\n# # ------------- Test Predictions\n# te_pred = np.clip(model.predict(X_te) + init_score, 0, 3)\n# te_pred_rounded = np.round(te_pred).astype(int)\n# for idx, pred in enumerate(te_pred_rounded):\n# te_preds[idx].append(pred)\n \n# # ------------- Feature Importance\n# if hasattr(model, \"feature_importance\"):\n# gains += model.feature_importance(importance_type='gain')\n# split += model.feature_importance(importance_type='split')\n \n# if hasattr(model, \"get_feature_importance\"):\n# feature_importances = model.get_feature_importance()\n# gains += np.array(feature_importances)\n# split += np.array(feature_importances) \n\n# trn_eval = cohen_kappa_score(y_tr, tr_pred_rounded, weights='quadratic')\n# val_eval = cohen_kappa_score(y_va, y_pred_rounded, weights='quadratic')\n# diff += abs(trn_eval - val_eval)\n\n# tr_preds_majority = np.array([mode(preds).mode if preds else 0 for preds in tr_preds])\n# oof_preds_majority = np.array([mode(preds).mode if preds else 0 for preds in oof_raw_preds])\n# te_preds_majority = np.array([mode(preds).mode if preds else 0 for preds in te_preds])\n\n# oof_score = cohen_kappa_score(y, oof_preds_majority, weights='quadratic')\n# tr_score = cohen_kappa_score(y, tr_preds_majority, weights='quadratic')\n\n# alpha = 0.8\n# combined_score = oof_score - alpha * (diff / n_splits)\n\n# return -combined_score","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.444649Z","iopub.execute_input":"2024-11-06T01:52:08.445039Z","iopub.status.idle":"2024-11-06T01:52:08.461985Z","shell.execute_reply.started":"2024-11-06T01:52:08.444999Z","shell.execute_reply":"2024-11-06T01:52:08.461003Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# study = optuna.create_study(direction=\"minimize\")\n# study.optimize(lambda trial: tune_cross_validate_model(trial, bst), n_trials=1000)\n\n# best_params = study.best_params\n# best_score = -study.best_value \n\n# print(\"Best parameters:\", best_params)\n# print(\"Best score:\", best_score)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.462971Z","iopub.execute_input":"2024-11-06T01:52:08.463315Z","iopub.status.idle":"2024-11-06T01:52:08.479418Z","shell.execute_reply.started":"2024-11-06T01:52:08.463280Z","shell.execute_reply":"2024-11-06T01:52:08.478273Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 100 runs\n# Best parameters: {'lambda_l2': 0.7676728933620794, 'feature_fraction': 0.49811372784537183, 'colsample_bytree': 0.10581494046323615, 'max_bin': 52, 'max_depth': 4, 'num_leaves': 19, 'min_data_in_leaf': 56, 'path_smooth': 2.5592238512022094}\n# Best score: 0.4804244004826107","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.480872Z","iopub.execute_input":"2024-11-06T01:52:08.481233Z","iopub.status.idle":"2024-11-06T01:52:08.490219Z","shell.execute_reply.started":"2024-11-06T01:52:08.481196Z","shell.execute_reply":"2024-11-06T01:52:08.489129Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"lgb_results['feature_importance']","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.491650Z","iopub.execute_input":"2024-11-06T01:52:08.492100Z","iopub.status.idle":"2024-11-06T01:52:08.507386Z","shell.execute_reply.started":"2024-11-06T01:52:08.492048Z","shell.execute_reply":"2024-11-06T01:52:08.506320Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def extended_cohen_kappa_score(labels, preds):\n f = np.sum((preds - labels)**2)\n g = np.sum((preds - a) ** 2 + b)\n return 1 - f / g ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.508576Z","iopub.execute_input":"2024-11-06T01:52:08.508929Z","iopub.status.idle":"2024-11-06T01:52:08.517131Z","shell.execute_reply.started":"2024-11-06T01:52:08.508890Z","shell.execute_reply":"2024-11-06T01:52:08.516082Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"oof_score = extended_cohen_kappa_score(\n y.to_numpy(),\n np.round(lgb_results['oof']).astype(int)\n)\n\ntr_score = extended_cohen_kappa_score(\n y.to_numpy(),\n np.round(lgb_results['trn']).astype(int)\n)\n\nprint(f\"{Fore.YELLOW}{Style.BRIGHT}# Extended TRAIN: tr_score={tr_score:.3f} {Style.RESET_ALL}\")\nprint(f\"{Fore.GREEN}{Style.BRIGHT}# Extended OOF: oof_score={oof_score:.3f} {Style.RESET_ALL}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.518731Z","iopub.execute_input":"2024-11-06T01:52:08.519374Z","iopub.status.idle":"2024-11-06T01:52:08.529676Z","shell.execute_reply.started":"2024-11-06T01:52:08.519305Z","shell.execute_reply":"2024-11-06T01:52:08.528636Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sub_test = test.to_pandas()\nsub_test['sii'] = np.round(lgb_results['test']).astype(int)\nsub_test['sii'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2024-11-06T01:52:08.531148Z","iopub.execute_input":"2024-11-06T01:52:08.531512Z","iopub.status.idle":"2024-11-06T01:52:08.549222Z","shell.execute_reply.started":"2024-11-06T01:52:08.531475Z","shell.execute_reply":"2024-11-06T01:52:08.548139Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sub_test[['id','sii']].to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2024-11-06T01:52:08.550513Z","iopub.execute_input":"2024-11-06T01:52:08.550942Z","iopub.status.idle":"2024-11-06T01:52:08.558098Z","shell.execute_reply.started":"2024-11-06T01:52:08.550892Z","shell.execute_reply":"2024-11-06T01:52:08.557074Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# TRAIN: tr_score=0.493 \n# OOF: oof_score=0.478 \n# LB = 0.455\n# OOF/LB = 0.466\n\n# TRAIN: tr_score=0.490 \n# OOF: oof_score=0.480 \n# LB = 0.449\n# OOF/LB = 0.464\n\n# TRAIN: tr_score=0.494 \n# OOF: oof_score=0.481 \n# LB = 0.462\n# OOF/LB = 0.0.472\n\n# TRAIN: tr_score=0.515 \n# OOF: oof_score=0.488 \n# LB = 0.463\n# OOF/LB = 0.476\n\n# TRAIN: tr_score=0.504 \n# OOF: oof_score=0.494 \n# LB = 0.457\n# OOF/LB = 0.476\n\n\n# TRAIN: 0.506\n# OOF: 0.498\n# DIFF: 0.038\n# LB = 0.461\n# OOF/LB = 0.480","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-06T01:52:08.559565Z","iopub.execute_input":"2024-11-06T01:52:08.559909Z","iopub.status.idle":"2024-11-06T01:52:08.569334Z","shell.execute_reply.started":"2024-11-06T01:52:08.559873Z","shell.execute_reply":"2024-11-06T01:52:08.568380Z"}},"outputs":[],"execution_count":null}]}