From 74ccafb7867992eaef445952ad63411023bfcc93 Mon Sep 17 00:00:00 2001 From: Minami Ueda Date: Fri, 27 Mar 2026 11:42:08 -0400 Subject: [PATCH 1/3] Setup ruff --- pyproject.toml | 24 ++++++++++++++++++++++++ uv.lock | 27 +++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 51ff646..59266ec 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,7 @@ dev = [ "pytest>=7.0.0", "pytest-cov>=4.0.0", "mypy>=1.0.0", + "ruff>=0.8.0", ] docs = [ "sphinx>=8.0.0; python_version >= '3.10'", @@ -80,3 +81,26 @@ exclude_lines = [ markers = [ "slow: marks tests as slow (deselect with '-m \"not slow\"')", ] + +[tool.ruff] +target-version = "py38" +line-length = 88 + +[tool.ruff.lint] +select = [ + "E", # pycodestyle errors + "W", # pycodestyle warnings + "F", # pyflakes + "I", # isort + "UP", # pyupgrade + "B", # flake8-bugbear +] +ignore = [ + "E501", # line length handled by ruff format +] + +[tool.ruff.lint.per-file-ignores] +"src/tailestim/estimators/plot/*" = ["F841", "E741"] + +[tool.ruff.lint.isort] +known-first-party = ["tailestim"] diff --git a/uv.lock b/uv.lock index a587c0d..246f7b9 100644 --- a/uv.lock +++ b/uv.lock @@ -2742,6 +2742,31 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/04/54/6f679c435d28e0a568d8e8a7c0a93a09010818634c3c3907fc98d8983770/roman_numerals-4.1.0-py3-none-any.whl", hash = "sha256:647ba99caddc2cc1e55a51e4360689115551bf4476d90e8162cf8c345fe233c7", size = 7676, upload-time = "2025-12-17T18:25:33.098Z" }, ] +[[package]] +name = "ruff" +version = "0.15.8" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/14/b0/73cf7550861e2b4824950b8b52eebdcc5adc792a00c514406556c5b80817/ruff-0.15.8.tar.gz", hash = "sha256:995f11f63597ee362130d1d5a327a87cb6f3f5eae3094c620bcc632329a4d26e", size = 4610921, upload-time = "2026-03-26T18:39:38.675Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/4a/92/c445b0cd6da6e7ae51e954939cb69f97e008dbe750cfca89b8cedc081be7/ruff-0.15.8-py3-none-linux_armv6l.whl", hash = "sha256:cbe05adeba76d58162762d6b239c9056f1a15a55bd4b346cfd21e26cd6ad7bc7", size = 10527394, upload-time = "2026-03-26T18:39:41.566Z" }, + { url = "https://files.pythonhosted.org/packages/eb/92/f1c662784d149ad1414cae450b082cf736430c12ca78367f20f5ed569d65/ruff-0.15.8-py3-none-macosx_10_12_x86_64.whl", hash = "sha256:d3e3d0b6ba8dca1b7ef9ab80a28e840a20070c4b62e56d675c24f366ef330570", size = 10905693, upload-time = "2026-03-26T18:39:30.364Z" }, + { url = "https://files.pythonhosted.org/packages/ca/f2/7a631a8af6d88bcef997eb1bf87cc3da158294c57044aafd3e17030613de/ruff-0.15.8-py3-none-macosx_11_0_arm64.whl", hash = "sha256:6ee3ae5c65a42f273f126686353f2e08ff29927b7b7e203b711514370d500de3", size = 10323044, upload-time = "2026-03-26T18:39:33.37Z" }, + { url = "https://files.pythonhosted.org/packages/67/18/1bf38e20914a05e72ef3b9569b1d5c70a7ef26cd188d69e9ca8ef588d5bf/ruff-0.15.8-py3-none-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:fdce027ada77baa448077ccc6ebb2fa9c3c62fd110d8659d601cf2f475858d94", size = 10629135, upload-time = "2026-03-26T18:39:44.142Z" }, + { url = "https://files.pythonhosted.org/packages/d2/e9/138c150ff9af60556121623d41aba18b7b57d95ac032e177b6a53789d279/ruff-0.15.8-py3-none-manylinux_2_17_armv7l.manylinux2014_armv7l.whl", hash = "sha256:12e617fc01a95e5821648a6df341d80456bd627bfab8a829f7cfc26a14a4b4a3", size = 10348041, upload-time = "2026-03-26T18:39:52.178Z" }, + { url = "https://files.pythonhosted.org/packages/02/f1/5bfb9298d9c323f842c5ddeb85f1f10ef51516ac7a34ba446c9347d898df/ruff-0.15.8-py3-none-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:432701303b26416d22ba696c39f2c6f12499b89093b61360abc34bcc9bf07762", size = 11121987, upload-time = "2026-03-26T18:39:55.195Z" }, + { url = "https://files.pythonhosted.org/packages/10/11/6da2e538704e753c04e8d86b1fc55712fdbdcc266af1a1ece7a51fff0d10/ruff-0.15.8-py3-none-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:d910ae974b7a06a33a057cb87d2a10792a3b2b3b35e33d2699fdf63ec8f6b17a", size = 11951057, upload-time = "2026-03-26T18:39:19.18Z" }, + { url = "https://files.pythonhosted.org/packages/83/f0/c9208c5fd5101bf87002fed774ff25a96eea313d305f1e5d5744698dc314/ruff-0.15.8-py3-none-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:2033f963c43949d51e6fdccd3946633c6b37c484f5f98c3035f49c27395a8ab8", size = 11464613, upload-time = "2026-03-26T18:40:06.301Z" }, + { url = "https://files.pythonhosted.org/packages/f8/22/d7f2fabdba4fae9f3b570e5605d5eb4500dcb7b770d3217dca4428484b17/ruff-0.15.8-py3-none-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:0f29b989a55572fb885b77464cf24af05500806ab4edf9a0fd8977f9759d85b1", size = 11257557, upload-time = "2026-03-26T18:39:57.972Z" }, + { url = "https://files.pythonhosted.org/packages/71/8c/382a9620038cf6906446b23ce8632ab8c0811b8f9d3e764f58bedd0c9a6f/ruff-0.15.8-py3-none-manylinux_2_31_riscv64.whl", hash = "sha256:ac51d486bf457cdc985a412fb1801b2dfd1bd8838372fc55de64b1510eff4bec", size = 11169440, upload-time = "2026-03-26T18:39:22.205Z" }, + { url = "https://files.pythonhosted.org/packages/4d/0d/0994c802a7eaaf99380085e4e40c845f8e32a562e20a38ec06174b52ef24/ruff-0.15.8-py3-none-musllinux_1_2_aarch64.whl", hash = "sha256:c9861eb959edab053c10ad62c278835ee69ca527b6dcd72b47d5c1e5648964f6", size = 10605963, upload-time = "2026-03-26T18:39:46.682Z" }, + { url = "https://files.pythonhosted.org/packages/19/aa/d624b86f5b0aad7cef6bbf9cd47a6a02dfdc4f72c92a337d724e39c9d14b/ruff-0.15.8-py3-none-musllinux_1_2_armv7l.whl", hash = "sha256:8d9a5b8ea13f26ae90838afc33f91b547e61b794865374f114f349e9036835fb", size = 10357484, upload-time = "2026-03-26T18:39:49.176Z" }, + { url = "https://files.pythonhosted.org/packages/35/c3/e0b7835d23001f7d999f3895c6b569927c4d39912286897f625736e1fd04/ruff-0.15.8-py3-none-musllinux_1_2_i686.whl", hash = "sha256:c2a33a529fb3cbc23a7124b5c6ff121e4d6228029cba374777bd7649cc8598b8", size = 10830426, upload-time = "2026-03-26T18:40:03.702Z" }, + { url = "https://files.pythonhosted.org/packages/f0/51/ab20b322f637b369383adc341d761eaaa0f0203d6b9a7421cd6e783d81b9/ruff-0.15.8-py3-none-musllinux_1_2_x86_64.whl", hash = "sha256:75e5cd06b1cf3f47a3996cfc999226b19aa92e7cce682dcd62f80d7035f98f49", size = 11345125, upload-time = "2026-03-26T18:39:27.799Z" }, + { url = "https://files.pythonhosted.org/packages/37/e6/90b2b33419f59d0f2c4c8a48a4b74b460709a557e8e0064cf33ad894f983/ruff-0.15.8-py3-none-win32.whl", hash = "sha256:bc1f0a51254ba21767bfa9a8b5013ca8149dcf38092e6a9eb704d876de94dc34", size = 10571959, upload-time = "2026-03-26T18:39:36.117Z" }, + { url = "https://files.pythonhosted.org/packages/1f/a2/ef467cb77099062317154c63f234b8a7baf7cb690b99af760c5b68b9ee7f/ruff-0.15.8-py3-none-win_amd64.whl", hash = "sha256:04f79eff02a72db209d47d665ba7ebcad609d8918a134f86cb13dd132159fc89", size = 11743893, upload-time = "2026-03-26T18:39:25.01Z" }, + { url = "https://files.pythonhosted.org/packages/15/e2/77be4fff062fa78d9b2a4dea85d14785dac5f1d0c1fb58ed52331f0ebe28/ruff-0.15.8-py3-none-win_arm64.whl", hash = "sha256:cf891fa8e3bb430c0e7fac93851a5978fc99c8fa2c053b57b118972866f8e5f2", size = 11048175, upload-time = "2026-03-26T18:40:01.06Z" }, +] + [[package]] name = "six" version = "1.17.0" @@ -2983,6 +3008,7 @@ dev = [ { name = "pytest", version = "9.0.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.10'" }, { name = "pytest-cov", version = "5.0.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.9'" }, { name = "pytest-cov", version = "7.1.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.9'" }, + { name = "ruff" }, ] docs = [ { name = "numpydoc", marker = "python_full_version >= '3.10'" }, @@ -3006,6 +3032,7 @@ dev = [ { name = "mypy", specifier = ">=1.0.0" }, { name = "pytest", specifier = ">=7.0.0" }, { name = "pytest-cov", specifier = ">=4.0.0" }, + { name = "ruff", specifier = ">=0.8.0" }, ] docs = [ { name = "numpydoc", marker = "python_full_version >= '3.10'", specifier = ">=1.5.0" }, From f4311af91070c8a9637e84352f7475d4c626ebd4 Mon Sep 17 00:00:00 2001 From: Minami Ueda Date: Fri, 27 Mar 2026 11:42:54 -0400 Subject: [PATCH 2/3] Add ruff linter/formatter and apply fixes --- src/tailestim/__init__.py | 24 +- src/tailestim/datasets.py | 44 +- src/tailestim/estimators/__init__.py | 18 +- src/tailestim/estimators/base.py | 57 +- src/tailestim/estimators/bulk_fit.py | 219 +-- src/tailestim/estimators/estimator_set.py | 125 +- src/tailestim/estimators/hill.py | 96 +- src/tailestim/estimators/kernel.py | 102 +- src/tailestim/estimators/moments.py | 100 +- src/tailestim/estimators/pickands.py | 40 +- src/tailestim/estimators/plot/plot_methods.py | 1208 +++++++++++------ src/tailestim/estimators/result.py | 53 +- src/tailestim/estimators/smooth_hill.py | 41 +- src/tailestim/estimators/tail_methods.py | 994 ++++++++------ tests/test_datasets.py | 71 +- tests/test_estimator_set.py | 97 +- tests/test_tail_methods.py | 277 +++- tests/test_tailestimation_validation.py | 7 +- 18 files changed, 2226 insertions(+), 1347 deletions(-) diff --git a/src/tailestim/__init__.py b/src/tailestim/__init__.py index ddb92a3..7c2c67c 100644 --- a/src/tailestim/__init__.py +++ b/src/tailestim/__init__.py @@ -3,23 +3,21 @@ # SPDX-License-Identifier: MIT from .datasets import TailData - from .estimators.base import BaseTailEstimator +from .estimators.estimator_set import TailEstimatorSet from .estimators.hill import HillEstimator -from .estimators.moments import MomentsEstimator from .estimators.kernel import KernelTypeEstimator +from .estimators.moments import MomentsEstimator from .estimators.pickands import PickandsEstimator from .estimators.smooth_hill import SmoothHillEstimator -from .estimators.estimator_set import TailEstimatorSet - __all__ = [ - 'TailData', - 'BaseTailEstimator', - 'HillEstimator', - 'MomentsEstimator', - 'KernelTypeEstimator', - 'PickandsEstimator', - 'SmoothHillEstimator', - 'TailEstimatorSet', -] \ No newline at end of file + "TailData", + "BaseTailEstimator", + "HillEstimator", + "MomentsEstimator", + "KernelTypeEstimator", + "PickandsEstimator", + "SmoothHillEstimator", + "TailEstimatorSet", +] diff --git a/src/tailestim/datasets.py b/src/tailestim/datasets.py index 556eead..bcb1a2c 100644 --- a/src/tailestim/datasets.py +++ b/src/tailestim/datasets.py @@ -1,13 +1,15 @@ +import logging import os + import numpy as np -import logging + class TailData: """Load and manage tail distribution datasets. - + This class provides functionality to load datasets either from the package's built-in data directory using a name, or from a custom path provided by the user. - + Parameters ---------- name : str, optional @@ -16,7 +18,7 @@ class TailData: path : str, optional Path to a custom dataset file. If provided, this takes precedence over `name`. Must be provided if `name` is None. - + Attributes ---------- name : str or None @@ -25,39 +27,39 @@ class TailData: Path to the dataset file if a custom dataset was loaded. data : numpy.ndarray The loaded dataset as a numpy array. - + Examples -------- Load a built-in dataset: - + >>> data = TailData(name='CAIDA_KONECT') >>> print(len(data.data)) - + Load a custom dataset: - + >>> data = TailData(path='path/to/my/data.dat') >>> print(len(data.data)) """ - + def __init__(self, name=None, path=None): if name is None and path is None: raise ValueError("Either 'name' or 'path' must be provided") - + if name is not None and path is not None: logging.info("Both 'name' and 'path' provided; 'path' will take precedence") - + self.name = name self.path = path self.data = self.load_data() def load_data(self): """Load data from either a built-in dataset or a custom file path. - + Returns ------- numpy.ndarray The loaded dataset as a numpy array. - + Raises ------ FileNotFoundError @@ -70,8 +72,8 @@ def load_data(self): logging.info(f"Using custom path: {file_path}") else: # Use the package data directory with the provided name - data_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data') - file_path = os.path.join(data_dir, f'{self.name}.dat') + data_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data") + file_path = os.path.join(data_dir, f"{self.name}.dat") logging.info(f"Using package data path: {file_path}") # Check if the file exists @@ -79,11 +81,13 @@ def load_data(self): if self.path is not None: raise FileNotFoundError(f"Data file not found at path: {file_path}") else: - raise FileNotFoundError(f"Data file '{self.name}.dat' not found in package data directory.") + raise FileNotFoundError( + f"Data file '{self.name}.dat' not found in package data directory." + ) # Load the data from the file using the provided method logging.info(f"Loading data from file: {file_path}") - with open(file_path, 'r') as file: + with open(file_path) as file: lines = file.readlines() # Determine the total number of data points @@ -94,14 +98,14 @@ def load_data(self): # Populate the ordered_data array for line in lines: degree, count = line.strip().split() - ordered_data[current_index:current_index + int(count)] = float(degree) + ordered_data[current_index : current_index + int(count)] = float(degree) current_index += int(count) return ordered_data def __repr__(self): """Return a string representation of the TailData object. - + Returns ------- str @@ -110,4 +114,4 @@ def __repr__(self): if self.path is not None: return f"TailData(path='{self.path}', data_length={len(self.data)})" else: - return f"TailData(name='{self.name}', data_length={len(self.data)})" \ No newline at end of file + return f"TailData(name='{self.name}', data_length={len(self.data)})" diff --git a/src/tailestim/estimators/__init__.py b/src/tailestim/estimators/__init__.py index 11f2bb2..a7023f9 100644 --- a/src/tailestim/estimators/__init__.py +++ b/src/tailestim/estimators/__init__.py @@ -1,19 +1,19 @@ """Estimators for tail index estimation.""" from .base import BaseTailEstimator +from .estimator_set import TailEstimatorSet from .hill import HillEstimator -from .moments import MomentsEstimator from .kernel import KernelTypeEstimator +from .moments import MomentsEstimator from .pickands import PickandsEstimator from .smooth_hill import SmoothHillEstimator -from .estimator_set import TailEstimatorSet __all__ = [ - 'BaseTailEstimator', - 'HillEstimator', - 'MomentsEstimator', - 'KernelTypeEstimator', - 'PickandsEstimator', - 'SmoothHillEstimator', - 'TailEstimatorSet', + "BaseTailEstimator", + "HillEstimator", + "MomentsEstimator", + "KernelTypeEstimator", + "PickandsEstimator", + "SmoothHillEstimator", + "TailEstimatorSet", ] diff --git a/src/tailestim/estimators/base.py b/src/tailestim/estimators/base.py index 126b228..f367d0b 100644 --- a/src/tailestim/estimators/base.py +++ b/src/tailestim/estimators/base.py @@ -1,10 +1,14 @@ """Base class for tail index estimation.""" -import numpy as np + from abc import ABC, abstractmethod -from typing import Dict, Any, Tuple, Union -from numpy.random import BitGenerator, SeedSequence, RandomState, Generator +from typing import Any, Dict, Tuple, Union + +import numpy as np +from numpy.random import BitGenerator, Generator, RandomState, SeedSequence + from .result import TailEstimatorResult + class BaseTailEstimator(ABC): """Abstract base class for tail index estimation. @@ -22,12 +26,15 @@ class BaseTailEstimator(ABC): **kwargs : dict Additional parameters specific to each estimation method. """ - + def __init__( - self, - bootstrap: bool = True, - base_seed: Union[None, SeedSequence, BitGenerator, Generator, RandomState] = None, - **kwargs): + self, + bootstrap: bool = True, + base_seed: Union[ + None, SeedSequence, BitGenerator, Generator, RandomState + ] = None, + **kwargs, + ): self.bootstrap = bootstrap self.base_seed = base_seed self.kwargs = kwargs @@ -36,12 +43,12 @@ def __init__( @abstractmethod def _estimate(self, ordered_data: np.ndarray) -> Tuple: """Core estimation method to be implemented by each specific estimator. - + Parameters ---------- ordered_data : np.ndarray Data array in decreasing order. - + Returns ------- Tuple @@ -51,34 +58,32 @@ def _estimate(self, ordered_data: np.ndarray) -> Tuple: def fit(self, data: np.ndarray) -> None: """Fit the estimator to the data. - + Parameters ---------- data : np.ndarray Input data array (e.g., degree sequence). The data will automatically be sorted in decreasing order. """ - ordered_data = np.sort(data)[::-1] # Each estimating functions require the data to be in decreasing order + ordered_data = np.sort(data)[ + ::-1 + ] # Each estimating functions require the data to be in decreasing order self.results = self._estimate(ordered_data) @abstractmethod def get_params(self) -> Dict[str, Any]: """Get the parameters of the estimator. - + Returns ------- dict Dictionary containing the parameters of the estimator. """ - return { - "bootstrap": self.bootstrap, - "base_seed": self.base_seed, - **self.kwargs - } + return {"bootstrap": self.bootstrap, "base_seed": self.base_seed, **self.kwargs} @abstractmethod def get_result(self) -> TailEstimatorResult: """Get the estimated parameters. - + Returns ------- TailEstimatorResult @@ -96,7 +101,7 @@ def get_result(self) -> TailEstimatorResult: if self.results is None: raise ValueError("Model not fitted yet. Call fit() first.") return TailEstimatorResult() - + def __repr__(self) -> str: """Return a string representation of the estimator.""" return f"{self.__class__.__name__}(bootstrap={self.bootstrap}, base_seed={self.base_seed}, kwargs={self.kwargs})" @@ -108,19 +113,19 @@ def __str__(self) -> str: estim_str += f"Estimator Type: {self.__class__.__name__}\n" estim_str += "-" * 50 + "\n" estim_str += f"Fitted: {'Yes' if self.results is not None else 'No'}\n" - + # Add the arguments provided estim_str += "Arguments:\n" estim_str += f" bootstrap: {self.bootstrap}\n" estim_str += f" base_seed: {self.base_seed}\n" - + # Add any additional kwargs if self.kwargs: for key, value in self.kwargs.items(): estim_str += f" {key}: {value}\n" - - # If the model is not fitted, return just the estim_str + + # If the model is not fitted, return just the estim_str if self.results is None: return estim_str + "Model not fitted yet. Call fit() first." - - return estim_str + + return estim_str diff --git a/src/tailestim/estimators/bulk_fit.py b/src/tailestim/estimators/bulk_fit.py index 06b003f..3868286 100644 --- a/src/tailestim/estimators/bulk_fit.py +++ b/src/tailestim/estimators/bulk_fit.py @@ -1,23 +1,37 @@ -import time -import numpy as np - import logging -logging.basicConfig(level=logging.WARNING) +import time -from .tail_methods import add_uniform_noise, get_distribution, get_ccdf +import numpy as np from .hill import HillEstimator -from .moments import MomentsEstimator from .kernel import KernelTypeEstimator +from .moments import MomentsEstimator from .pickands import PickandsEstimator from .smooth_hill import SmoothHillEstimator +from .tail_methods import add_uniform_noise, get_ccdf, get_distribution + +logging.basicConfig(level=logging.WARNING) + -def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hsteps=200, - bootstrap_flag=True, t_bootstrap=0.5, r_bootstrap=500, diagn_plots=False, - eps_stop=1.0, noise_flag=True, verbose=False, p_noise=1, base_seed=None): +def fit_estimators( + ordered_data, + number_of_bins=30, + r_smooth=2, + alpha=0.6, + hsteps=200, + bootstrap_flag=True, + t_bootstrap=0.5, + r_bootstrap=500, + diagn_plots=False, + eps_stop=1.0, + noise_flag=True, + verbose=False, + p_noise=1, + base_seed=None, +): """ Fit various tail estimators to the data at once. - + Parameters ---------- ordered_data : numpy.ndarray @@ -39,7 +53,7 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep diagn_plots : bool, default=False Switch on/off generation of AMSE diagnostic plots. eps_stop : float, default=1.0 - Parameter controlling range of AMSE minimization. Defined as the fraction of order + Parameter controlling range of AMSE minimization. Defined as the fraction of order statistics to consider during the AMSE minimization step. noise_flag : bool, default=True Switch on/off uniform noise in range [-5*10^(-p), 5*10^(-p)] that is added to each @@ -50,7 +64,7 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep Parameter controlling noise amplitude. base_seed : int, optional Base random seed for reproducibility of bootstrap. Only used for methods with bootstrap. - + Returns ------- dict @@ -58,33 +72,32 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep different estimators (Pickands, Hill, smooth Hill, moments, kernel-type). """ results = {} - + # calculate log-binned PDF logging.debug("Calculating PDF...") t1 = time.time() x_pdf, y_pdf = get_distribution(ordered_data, number_of_bins=number_of_bins) t2 = time.time() - logging.debug("Elapsed time(PDF):", t2-t1) - results['pdf'] = {'x': x_pdf, 'y': y_pdf} + logging.debug("Elapsed time(PDF):", t2 - t1) + results["pdf"] = {"x": x_pdf, "y": y_pdf} # calculate CCDF logging.debug("Calculating CCDF...") t1 = time.time() x_ccdf, y_ccdf = get_ccdf(ordered_data) t2 = time.time() - logging.debug("Elapsed time:", t2-t1) - results['ccdf'] = {'x': x_ccdf, 'y': y_ccdf} + logging.debug("Elapsed time:", t2 - t1) + results["ccdf"] = {"x": x_ccdf, "y": y_ccdf} # add noise if needed if noise_flag: - original_discrete_data = ordered_data discrete_ordered_data = ordered_data.copy() discrete_ordered_data[::-1].sort() ordered_data = add_uniform_noise(ordered_data, p=p_noise, base_seed=base_seed) ordered_data[::-1].sort() - results['ordered_data'] = ordered_data + results["ordered_data"] = ordered_data if noise_flag: - results['discrete_ordered_data'] = discrete_ordered_data + results["discrete_ordered_data"] = discrete_ordered_data # perform Pickands estimation logging.debug("Calculating Pickands...") @@ -94,8 +107,8 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep pickands_result = pickands.get_result() k_p_arr, xi_p_arr = pickands_result.k_arr_, pickands_result.xi_arr_ t2 = time.time() - logging.debug("Elapsed time (Pickands):", t2-t1) - results['pickands'] = {'k_arr_': k_p_arr, 'xi_arr_': xi_p_arr} + logging.debug("Elapsed time (Pickands):", t2 - t1) + results["pickands"] = {"k_arr_": k_p_arr, "xi_arr_": xi_p_arr} # perform smooth Hill estimation logging.debug("Calculating smooth Hill...") @@ -105,8 +118,8 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep smooth_hill_result = smooth_hill.get_result() k_sh_arr, xi_sh_arr = smooth_hill_result.k_arr_, smooth_hill_result.xi_arr_ t2 = time.time() - logging.debug("Elapsed time (smooth Hill):", t2-t1) - results['smooth_hill'] = {'k_arr_': k_sh_arr, 'xi_arr_': xi_sh_arr} + logging.debug("Elapsed time (smooth Hill):", t2 - t1) + results["smooth_hill"] = {"k_arr_": k_sh_arr, "xi_arr_": xi_sh_arr} # perform adjusted Hill estimation logging.debug("Calculating adjusted Hill...") @@ -118,16 +131,16 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep diagn_plots=diagn_plots, eps_stop=eps_stop, verbose=verbose, - base_seed=base_seed + base_seed=base_seed, ) hill.fit(ordered_data) hill_result = hill.get_result() k_h_arr = hill_result.k_arr_ xi_h_arr = hill_result.xi_arr_ - - results['hill'] = {'k_arr_': k_h_arr, 'xi_arr_': xi_h_arr} - - if bootstrap_flag and hasattr(hill_result, 'k_star_'): + + results["hill"] = {"k_arr_": k_h_arr, "xi_arr_": xi_h_arr} + + if bootstrap_flag and hasattr(hill_result, "k_star_"): k_h_star = hill_result.k_star_ xi_h_star = hill_result.xi_star_ x1_h_arr = hill_result.bootstrap_results_.first_bootstrap_.x_arr_ @@ -138,29 +151,31 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep n2_h_amse = hill_result.bootstrap_results_.second_bootstrap_.amse_ k2_h = hill_result.bootstrap_results_.second_bootstrap_.k_min_ max_h_index2 = hill_result.bootstrap_results_.second_bootstrap_.max_index_ - - results['hill'].update({ - 'k_star_': k_h_star, - 'xi_star_': xi_h_star, - 'bootstrap_results_': { - 'first_bootstrap_': { - 'x_arr_': x1_h_arr, - 'amse_': n1_h_amse, - 'k_min_': k1_h, - 'max_index_': max_h_index1 + + results["hill"].update( + { + "k_star_": k_h_star, + "xi_star_": xi_h_star, + "bootstrap_results_": { + "first_bootstrap_": { + "x_arr_": x1_h_arr, + "amse_": n1_h_amse, + "k_min_": k1_h, + "max_index_": max_h_index1, + }, + "second_bootstrap_": { + "x_arr_": x2_h_arr, + "amse_": n2_h_amse, + "k_min_": k2_h, + "max_index_": max_h_index2, + }, }, - 'second_bootstrap_': { - 'x_arr_': x2_h_arr, - 'amse_': n2_h_amse, - 'k_min_': k2_h, - 'max_index_': max_h_index2 - } } - }) - + ) + t2 = time.time() if verbose: - logging.debug("Elapsed time (Hill):", t2-t1) + logging.debug("Elapsed time (Hill):", t2 - t1) # perform moments estimation if verbose: @@ -173,16 +188,16 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep diagn_plots=diagn_plots, eps_stop=eps_stop, verbose=verbose, - base_seed=base_seed + base_seed=base_seed, ) moments.fit(ordered_data) moments_result = moments.get_result() k_m_arr = moments_result.k_arr_ xi_m_arr = moments_result.xi_arr_ - - results['moments'] = {'k_arr_': k_m_arr, 'xi_arr_': xi_m_arr} - - if bootstrap_flag and hasattr(moments_result, 'k_star_'): + + results["moments"] = {"k_arr_": k_m_arr, "xi_arr_": xi_m_arr} + + if bootstrap_flag and hasattr(moments_result, "k_star_"): k_m_star = moments_result.k_star_ xi_m_star = moments_result.xi_star_ x1_m_arr = moments_result.bootstrap_results_.first_bootstrap_.x_arr_ @@ -193,29 +208,31 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep n2_m_amse = moments_result.bootstrap_results_.second_bootstrap_.amse_ k2_m = moments_result.bootstrap_results_.second_bootstrap_.k_min_ max_m_index2 = moments_result.bootstrap_results_.second_bootstrap_.max_index_ - - results['moments'].update({ - 'k_star_': k_m_star, - 'xi_star_': xi_m_star, - 'bootstrap_results_': { - 'first_bootstrap_': { - 'x_arr_': x1_m_arr, - 'amse_': n1_m_amse, - 'k_min_': k1_m, - 'max_index_': max_m_index1 + + results["moments"].update( + { + "k_star_": k_m_star, + "xi_star_": xi_m_star, + "bootstrap_results_": { + "first_bootstrap_": { + "x_arr_": x1_m_arr, + "amse_": n1_m_amse, + "k_min_": k1_m, + "max_index_": max_m_index1, + }, + "second_bootstrap_": { + "x_arr_": x2_m_arr, + "amse_": n2_m_amse, + "k_min_": k2_m, + "max_index_": max_m_index2, + }, }, - 'second_bootstrap_': { - 'x_arr_': x2_m_arr, - 'amse_': n2_m_amse, - 'k_min_': k2_m, - 'max_index_': max_m_index2 - } } - }) - + ) + t2 = time.time() if verbose: - logging.debug("Elapsed time (moments):", t2-t1) + logging.debug("Elapsed time (moments):", t2 - t1) # perform kernel-type estimation if verbose: @@ -230,16 +247,16 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep diagn_plots=diagn_plots, eps_stop=eps_stop, verbose=verbose, - base_seed=base_seed + base_seed=base_seed, ) kernel.fit(ordered_data) kernel_result = kernel.get_result() k_k_arr = kernel_result.k_arr_ xi_k_arr = kernel_result.xi_arr_ - - results['kernel'] = {'k_arr_': k_k_arr, 'xi_arr_': xi_k_arr} - - if bootstrap_flag and hasattr(kernel_result, 'k_star_'): + + results["kernel"] = {"k_arr_": k_k_arr, "xi_arr_": xi_k_arr} + + if bootstrap_flag and hasattr(kernel_result, "k_star_"): k_k_star = kernel_result.k_star_ xi_k_star = kernel_result.xi_star_ x1_k_arr = kernel_result.bootstrap_results_.first_bootstrap_.x_arr_ @@ -250,33 +267,35 @@ def fit_estimators(ordered_data, number_of_bins=30, r_smooth=2, alpha=0.6, hstep n2_k_amse = kernel_result.bootstrap_results_.second_bootstrap_.amse_ h2 = kernel_result.bootstrap_results_.second_bootstrap_.h_min_ max_k_index2 = kernel_result.bootstrap_results_.second_bootstrap_.max_index_ - - results['kernel'].update({ - 'k_star_': k_k_star, - 'xi_star_': xi_k_star, - 'bootstrap_results_': { - 'first_bootstrap_': { - 'x_arr_': x1_k_arr, - 'amse_': n1_k_amse, - 'h_min': h1, - 'max_index_': max_k_index1 + + results["kernel"].update( + { + "k_star_": k_k_star, + "xi_star_": xi_k_star, + "bootstrap_results_": { + "first_bootstrap_": { + "x_arr_": x1_k_arr, + "amse_": n1_k_amse, + "h_min": h1, + "max_index_": max_k_index1, + }, + "second_bootstrap_": { + "x_arr_": x2_k_arr, + "amse_": n2_k_amse, + "h_min": h2, + "max_index_": max_k_index2, + }, }, - 'second_bootstrap_': { - 'x_arr_': x2_k_arr, - 'amse_': n2_k_amse, - 'h_min': h2, - 'max_index_': max_k_index2 - } } - }) - + ) + # Find the index of k_k_star in k_k_arr for plotting if k_k_star is not None: k_k1_star = np.argmin(np.abs(k_k_arr - k_k_star)) - results['kernel']['k_k1_star'] = k_k1_star - + results["kernel"]["k_k1_star"] = k_k1_star + t2 = time.time() if verbose: - logging.debug("Elapsed time (kernel-type):", t2-t1) - - return results \ No newline at end of file + logging.debug("Elapsed time (kernel-type):", t2 - t1) + + return results diff --git a/src/tailestim/estimators/estimator_set.py b/src/tailestim/estimators/estimator_set.py index 5bf0add..612f632 100644 --- a/src/tailestim/estimators/estimator_set.py +++ b/src/tailestim/estimators/estimator_set.py @@ -1,14 +1,17 @@ +from typing import Any, Dict, Optional, Tuple, Union + import numpy as np from matplotlib import pyplot as plt -from typing import Optional, Tuple, Union, Dict, Any, List -from numpy.random import BitGenerator, SeedSequence, RandomState, Generator +from numpy.random import BitGenerator, Generator, RandomState, SeedSequence + from .bulk_fit import fit_estimators from .plot.plot_methods import make_plots + class TailEstimatorSet: """ Class for running estimation with multiple estimator methods at once and creating a plot for comparison. - + Parameters ---------- data : np.ndarray @@ -57,10 +60,10 @@ class TailEstimatorSet: Whether to create the plots immediately upon initialization. base_seed: None | SeedSequence | BitGenerator | Generator | RandomState, default=None Base random seed for reproducibility of bootstrap. Only used for methods with bootstrap. - + """ - + def __init__( self, data: np.ndarray = None, @@ -81,7 +84,9 @@ def __init__( p_noise: int = 1, savedata: bool = False, auto_plot: bool = False, - base_seed: Union[None, SeedSequence, BitGenerator, Generator, RandomState] = None + base_seed: Union[ + None, SeedSequence, BitGenerator, Generator, RandomState + ] = None, ): # Store parameters self.output_file_path = output_file_path @@ -101,32 +106,32 @@ def __init__( self.p_noise = p_noise self.savedata = savedata self.base_seed = base_seed - + # Initialize data-related attributes self.data = None self.ordered_data = None self.results = None - + # Store figure and axes as None initially self.fig = None self.axes = None - + # Fit data if provided if data is not None: self.fit(data) - + # Create the plots immediately if auto_plot is True and data is provided if auto_plot and data is not None: self.plot() - - def fit(self, data: np.ndarray) -> 'TailEstimatorSet': + + def fit(self, data: np.ndarray) -> "TailEstimatorSet": """Fit the estimators to the data. - + Parameters ---------- data : np.ndarray The data to fit the estimators to. - + Returns ------- self : TailEstimatorSet @@ -134,11 +139,11 @@ def fit(self, data: np.ndarray) -> 'TailEstimatorSet': """ # Make sure data is a numpy array data_array = np.asarray(data) - + # Store the data self.data = data_array self.ordered_data = np.sort(data_array)[::-1] - + # Fit the estimators self.results = fit_estimators( ordered_data=self.ordered_data, @@ -154,18 +159,18 @@ def fit(self, data: np.ndarray) -> 'TailEstimatorSet': verbose=self.verbose, noise_flag=self.noise_flag, p_noise=self.p_noise, - base_seed=self.base_seed + base_seed=self.base_seed, ) - + # Reset figure and axes self.fig = None self.axes = None - + return self - + def plot(self) -> Tuple[plt.Figure, np.ndarray]: """Create and return the plots. - + Returns ------- fig : matplotlib.figure.Figure @@ -175,20 +180,20 @@ def plot(self) -> Tuple[plt.Figure, np.ndarray]: """ if self.ordered_data is None: raise ValueError("No data has been fitted. Call fit() first.") - + self.fig, self.axes = self._create_plots() return self.fig, self.axes - + def plot_diagnostics(self) -> Tuple[plt.Figure, np.ndarray]: """Create and return the diagnostic plots. - + Returns ------- fig_d : matplotlib.figure.Figure The diagnostic figure object. axes_d : numpy.ndarray Array of diagnostic axes objects. - + Raises ------ ValueError @@ -196,18 +201,22 @@ def plot_diagnostics(self) -> Tuple[plt.Figure, np.ndarray]: """ if self.ordered_data is None: raise ValueError("No data has been fitted. Call fit() first.") - + if not self.bootstrap_flag: - raise ValueError("Diagnostic plots require bootstrap to be enabled. Set bootstrap_flag=True when creating the TailEstimatorSet.") - + raise ValueError( + "Diagnostic plots require bootstrap to be enabled. Set bootstrap_flag=True when creating the TailEstimatorSet." + ) + if not self.diagnostic_plots: - raise ValueError("Diagnostic plots are not enabled. Set diagnostic_plots=True when creating the TailEstimatorSet.") - + raise ValueError( + "Diagnostic plots are not enabled. Set diagnostic_plots=True when creating the TailEstimatorSet." + ) + return self._create_diagnostic_plots() - + def _create_diagnostic_plots(self) -> Tuple[plt.Figure, np.ndarray]: """Create the diagnostic plots using the make_diagnostic_plots function. - + Returns ------- fig_d : matplotlib.figure.Figure @@ -216,7 +225,7 @@ def _create_diagnostic_plots(self) -> Tuple[plt.Figure, np.ndarray]: Array of diagnostic axes objects. """ from .plot.plot_methods import make_diagnostic_plots - + return make_diagnostic_plots( results=self.results, output_file_path=self.output_file_path, @@ -226,10 +235,10 @@ def _create_diagnostic_plots(self) -> Tuple[plt.Figure, np.ndarray]: noise_flag=self.noise_flag, savedata=self.savedata, ) - + def _create_plots(self) -> Tuple[plt.Figure, np.ndarray]: """Create the plots using the make_plots function. - + Returns ------- fig : matplotlib.figure.Figure @@ -250,47 +259,47 @@ def _create_plots(self) -> Tuple[plt.Figure, np.ndarray]: noise_flag=self.noise_flag, savedata=self.savedata, ) - + def get_params(self) -> Dict[str, Any]: """Get the parameters used for plotting. - + Returns ------- Dict[str, Any] Dictionary of parameters used for plotting. """ return { - 'data_length': len(self.data) if self.data is not None else 0, - 'number_of_bins': self.number_of_bins, - 'r_smooth': self.r_smooth, - 'alpha': self.alpha, - 'hsteps': self.hsteps, - 'bootstrap_flag': self.bootstrap_flag, - 't_bootstrap': self.t_bootstrap, - 'r_bootstrap': self.r_bootstrap, - 'diagnostic_plots': self.diagnostic_plots, - 'eps_stop': self.eps_stop, - 'theta1': self.theta1, - 'theta2': self.theta2, - 'verbose': self.verbose, - 'noise_flag': self.noise_flag, - 'p_noise': self.p_noise, - 'savedata': self.savedata, - 'base_seed': self.base_seed + "data_length": len(self.data) if self.data is not None else 0, + "number_of_bins": self.number_of_bins, + "r_smooth": self.r_smooth, + "alpha": self.alpha, + "hsteps": self.hsteps, + "bootstrap_flag": self.bootstrap_flag, + "t_bootstrap": self.t_bootstrap, + "r_bootstrap": self.r_bootstrap, + "diagnostic_plots": self.diagnostic_plots, + "eps_stop": self.eps_stop, + "theta1": self.theta1, + "theta2": self.theta2, + "verbose": self.verbose, + "noise_flag": self.noise_flag, + "p_noise": self.p_noise, + "savedata": self.savedata, + "base_seed": self.base_seed, } - + def __repr__(self) -> str: """Return a string representation of the object.""" if self.data is None: return "TailEstimatorSet(not fitted)" else: return f"TailEstimatorSet(data_length={len(self.data)})" - + def __call__(self) -> Tuple[plt.Figure, np.ndarray]: """Return the figure and axes when the object is called.""" if self.ordered_data is None: raise ValueError("No data has been fitted. Call fit() first.") - + if self.fig is None: self.plot() - return self.fig, self.axes \ No newline at end of file + return self.fig, self.axes diff --git a/src/tailestim/estimators/hill.py b/src/tailestim/estimators/hill.py index fc27e98..b737977 100644 --- a/src/tailestim/estimators/hill.py +++ b/src/tailestim/estimators/hill.py @@ -1,17 +1,21 @@ """Hill estimator implementation for tail index estimation.""" + +from typing import Any, Dict, Tuple, Union + import numpy as np -from typing import Dict, Any, Tuple, Union -from numpy.random import BitGenerator, SeedSequence, RandomState, Generator +from numpy.random import BitGenerator, Generator, RandomState, SeedSequence + from .base import BaseTailEstimator from .result import TailEstimatorResult from .tail_methods import hill_estimator as hill_estimate + class HillEstimator(BaseTailEstimator): """Hill estimator for tail index estimation. - + This class implements the Hill estimator with optional double-bootstrap for optimal threshold selection. - + Parameters ---------- bootstrap : bool, default=True @@ -45,9 +49,11 @@ def __init__( eps_stop: float = 0.99, verbose: bool = False, diagn_plots: bool = False, - base_seed: Union[None, SeedSequence, BitGenerator, Generator, RandomState] = None, + base_seed: Union[ + None, SeedSequence, BitGenerator, Generator, RandomState + ] = None, max_resample: int = 50, - **kwargs + **kwargs, ): super().__init__(bootstrap=bootstrap, base_seed=base_seed, **kwargs) self.t_bootstrap = t_bootstrap @@ -59,12 +65,12 @@ def __init__( def _estimate(self, ordered_data: np.ndarray) -> Tuple: """Estimate the tail index using the Hill estimator. - + Parameters ---------- ordered_data : np.ndarray Data array in decreasing order. - + Returns ------- Tuple @@ -79,12 +85,12 @@ def _estimate(self, ordered_data: np.ndarray) -> Tuple: diagn_plots=self.diagn_plots, eps_stop=self.eps_stop, base_seed=self.base_seed, - max_resample=self.max_resample + max_resample=self.max_resample, ) def get_params(self) -> Dict[str, Any]: """Get the parameters of the estimator. - + Returns ------- dict @@ -99,7 +105,7 @@ def get_params(self) -> Dict[str, Any]: "verbose": self.verbose, "diagn_plots": self.diagn_plots, "max_resample": self.max_resample, - **self.kwargs + **self.kwargs, } def get_result(self) -> TailEstimatorResult: @@ -129,36 +135,50 @@ def get_result(self) -> TailEstimatorResult: if self.results is None: raise ValueError("Model not fitted yet. Call fit() first.") - - k_arr, xi_arr, k_star, xi_star, x1_arr, n1_amse, k1, max_index1, \ - x2_arr, n2_amse, k2, max_index2 = self.results - + + ( + k_arr, + xi_arr, + k_star, + xi_star, + x1_arr, + n1_amse, + k1, + max_index1, + x2_arr, + n2_amse, + k2, + max_index2, + ) = self.results + res = { - 'k_arr_': k_arr, - 'xi_arr_': xi_arr, + "k_arr_": k_arr, + "xi_arr_": xi_arr, } - + if self.bootstrap and k_star is not None: - gamma = 1 + 1./xi_star - res.update({ - 'estimator': self, - 'k_star_': k_star, - 'xi_star_': xi_star, - 'gamma_': gamma, - 'bootstrap_results_': { - 'first_bootstrap_': { - 'x_arr_': x1_arr, - 'amse_': n1_amse, - 'k_min_': k1, - 'max_index_': max_index1 + gamma = 1 + 1.0 / xi_star + res.update( + { + "estimator": self, + "k_star_": k_star, + "xi_star_": xi_star, + "gamma_": gamma, + "bootstrap_results_": { + "first_bootstrap_": { + "x_arr_": x1_arr, + "amse_": n1_amse, + "k_min_": k1, + "max_index_": max_index1, + }, + "second_bootstrap_": { + "x_arr_": x2_arr, + "amse_": n2_amse, + "k_min_": k2, + "max_index_": max_index2, + }, }, - 'second_bootstrap_': { - 'x_arr_': x2_arr, - 'amse_': n2_amse, - 'k_min_': k2, - 'max_index_': max_index2 - } } - }) - + ) + return TailEstimatorResult(res) diff --git a/src/tailestim/estimators/kernel.py b/src/tailestim/estimators/kernel.py index a75b0db..cdf1cab 100644 --- a/src/tailestim/estimators/kernel.py +++ b/src/tailestim/estimators/kernel.py @@ -1,18 +1,22 @@ """Kernel-type estimator implementation for tail index estimation.""" + +from typing import Any, Dict, Tuple, Union + import numpy as np -from typing import Dict, Any, Tuple, Union -from numpy.random import BitGenerator, SeedSequence, RandomState, Generator +from numpy.random import BitGenerator, Generator, RandomState, SeedSequence + from .base import BaseTailEstimator from .result import TailEstimatorResult from .tail_methods import kernel_type_estimator as kernel_estimate + class KernelTypeEstimator(BaseTailEstimator): """Kernel-type estimator for tail index estimation. - + This class implements the Kernel-type estimator with optional double-bootstrap for optimal bandwidth selection. It uses both biweight and triweight kernels for estimation. - + Parameters ---------- bootstrap : bool, default=True @@ -38,7 +42,7 @@ class KernelTypeEstimator(BaseTailEstimator): base_seed: None | SeedSequence | BitGenerator | Generator | RandomState, default=None Base random seed for reproducibility of bootstrap. """ - + def __init__( self, bootstrap: bool = True, @@ -49,8 +53,10 @@ def __init__( eps_stop: float = 0.99, verbose: bool = False, diagn_plots: bool = False, - base_seed: Union[None, SeedSequence, BitGenerator, Generator, RandomState] = None, - **kwargs + base_seed: Union[ + None, SeedSequence, BitGenerator, Generator, RandomState + ] = None, + **kwargs, ): super().__init__(bootstrap=bootstrap, base_seed=base_seed, **kwargs) self.hsteps = hsteps @@ -63,12 +69,12 @@ def __init__( def _estimate(self, ordered_data: np.ndarray) -> Tuple: """Estimate tail index using kernel-type estimator. - + Parameters ---------- ordered_data : np.ndarray Data array in decreasing order. - + Returns ------- Tuple @@ -84,12 +90,12 @@ def _estimate(self, ordered_data: np.ndarray) -> Tuple: verbose=self.verbose, diagn_plots=self.diagn_plots, eps_stop=self.eps_stop, - base_seed=self.base_seed + base_seed=self.base_seed, ) def get_params(self) -> Dict[str, Any]: """Get the parameters of the estimator. - + Returns ------- dict @@ -105,7 +111,7 @@ def get_params(self) -> Dict[str, Any]: "diagn_plots": self.diagn_plots, "alpha": self.alpha, "hsteps": self.hsteps, - **self.kwargs + **self.kwargs, } def get_result(self) -> TailEstimatorResult: @@ -127,43 +133,57 @@ def get_result(self) -> TailEstimatorResult: Optimal order statistic (k*). bootstrap_results_ : dict Bootstrap results. - + Returns ------- TailEstimatorResult """ if self.results is None: raise ValueError("Model not fitted yet. Call fit() first.") - - k_arr, xi_arr, k_star, xi_star, x1_arr, n1_amse, h1, max_index1, \ - x2_arr, n2_amse, h2, max_index2 = self.results - + + ( + k_arr, + xi_arr, + k_star, + xi_star, + x1_arr, + n1_amse, + h1, + max_index1, + x2_arr, + n2_amse, + h2, + max_index2, + ) = self.results + res = { - 'k_arr_': k_arr, - 'xi_arr_': xi_arr, + "k_arr_": k_arr, + "xi_arr_": xi_arr, } - + if self.bootstrap and k_star is not None: - gamma = float('inf') if xi_star <= 0 else 1 + 1./xi_star - res.update({ - 'estimator': self, - 'k_star_': k_star, - 'xi_star_': xi_star, - 'gamma_': gamma, - 'bootstrap_results_': { - 'first_bootstrap_': { - 'x_arr_': x1_arr, - 'amse_': n1_amse, - 'h_min_': h1, - 'max_index_': max_index1 + gamma = float("inf") if xi_star <= 0 else 1 + 1.0 / xi_star + res.update( + { + "estimator": self, + "k_star_": k_star, + "xi_star_": xi_star, + "gamma_": gamma, + "bootstrap_results_": { + "first_bootstrap_": { + "x_arr_": x1_arr, + "amse_": n1_amse, + "h_min_": h1, + "max_index_": max_index1, + }, + "second_bootstrap_": { + "x_arr_": x2_arr, + "amse_": n2_amse, + "h_min_": h2, + "max_index_": max_index2, + }, }, - 'second_bootstrap_': { - 'x_arr_': x2_arr, - 'amse_': n2_amse, - 'h_min_': h2, - 'max_index_': max_index2 - } } - }) - - return TailEstimatorResult(res) \ No newline at end of file + ) + + return TailEstimatorResult(res) diff --git a/src/tailestim/estimators/moments.py b/src/tailestim/estimators/moments.py index 776f04e..4cde7fb 100644 --- a/src/tailestim/estimators/moments.py +++ b/src/tailestim/estimators/moments.py @@ -1,17 +1,21 @@ """Moments estimator implementation for tail index estimation.""" + +from typing import Any, Dict, Tuple, Union + import numpy as np -from typing import Dict, Any, Tuple, Union -from .result import TailEstimatorResult -from numpy.random import BitGenerator, SeedSequence, RandomState, Generator +from numpy.random import BitGenerator, Generator, RandomState, SeedSequence + from .base import BaseTailEstimator +from .result import TailEstimatorResult from .tail_methods import moments_estimator as moments_estimate + class MomentsEstimator(BaseTailEstimator): """Moments estimator for tail index estimation. - + This class implements the Moments estimator with optional double-bootstrap for optimal threshold selection. - + Parameters ---------- bootstrap : bool, default=True @@ -32,7 +36,7 @@ class MomentsEstimator(BaseTailEstimator): base_seed: None | SeedSequence | BitGenerator | Generator | RandomState, default=None Base random seed for reproducibility of bootstrap. """ - + def __init__( self, bootstrap: bool = True, @@ -41,8 +45,10 @@ def __init__( eps_stop: float = 0.99, verbose: bool = False, diagn_plots: bool = False, - base_seed: Union[None, SeedSequence, BitGenerator, Generator, RandomState] = None, - **kwargs + base_seed: Union[ + None, SeedSequence, BitGenerator, Generator, RandomState + ] = None, + **kwargs, ): super().__init__(bootstrap=bootstrap, base_seed=base_seed, **kwargs) self.t_bootstrap = t_bootstrap @@ -53,12 +59,12 @@ def __init__( def _estimate(self, ordered_data: np.ndarray) -> Tuple: """Estimate tail index using the Moments method. - + Parameters ---------- ordered_data : np.ndarray Data array in decreasing order. - + Returns ------- Tuple @@ -72,12 +78,12 @@ def _estimate(self, ordered_data: np.ndarray) -> Tuple: verbose=self.verbose, diagn_plots=self.diagn_plots, eps_stop=self.eps_stop, - base_seed=self.base_seed + base_seed=self.base_seed, ) def get_params(self) -> Dict[str, Any]: """Get the parameters of the estimator. - + Returns ------- dict @@ -91,7 +97,7 @@ def get_params(self) -> Dict[str, Any]: "eps_stop": self.eps_stop, "verbose": self.verbose, "diagn_plots": self.diagn_plots, - **self.kwargs + **self.kwargs, } def get_result(self) -> TailEstimatorResult: @@ -120,35 +126,49 @@ def get_result(self) -> TailEstimatorResult: """ if self.results is None: raise ValueError("Model not fitted yet. Call fit() first.") - - k_arr, xi_arr, k_star, xi_star, x1_arr, n1_amse, k1, max_index1, \ - x2_arr, n2_amse, k2, max_index2 = self.results - + + ( + k_arr, + xi_arr, + k_star, + xi_star, + x1_arr, + n1_amse, + k1, + max_index1, + x2_arr, + n2_amse, + k2, + max_index2, + ) = self.results + res = { - 'k_arr_': k_arr, - 'xi_arr_': xi_arr, + "k_arr_": k_arr, + "xi_arr_": xi_arr, } - + if self.bootstrap and k_star is not None: - gamma = float('inf') if xi_star <= 0 else 1 + 1./xi_star - res.update({ - 'estimator': self, - 'k_star_': k_star, - 'xi_star_': xi_star, - 'gamma_': gamma, - 'bootstrap_results_': { - 'first_bootstrap_': { - 'x_arr_': x1_arr, - 'amse_': n1_amse, - 'k_min_': k1, - 'max_index_': max_index1 + gamma = float("inf") if xi_star <= 0 else 1 + 1.0 / xi_star + res.update( + { + "estimator": self, + "k_star_": k_star, + "xi_star_": xi_star, + "gamma_": gamma, + "bootstrap_results_": { + "first_bootstrap_": { + "x_arr_": x1_arr, + "amse_": n1_amse, + "k_min_": k1, + "max_index_": max_index1, + }, + "second_bootstrap_": { + "x_arr_": x2_arr, + "amse_": n2_amse, + "k_min_": k2, + "max_index_": max_index2, + }, }, - 'second_bootstrap_': { - 'x_arr_': x2_arr, - 'amse_': n2_amse, - 'k_min_': k2, - 'max_index_': max_index2 - } } - }) - return TailEstimatorResult(res) \ No newline at end of file + ) + return TailEstimatorResult(res) diff --git a/src/tailestim/estimators/pickands.py b/src/tailestim/estimators/pickands.py index 090ba73..141596b 100644 --- a/src/tailestim/estimators/pickands.py +++ b/src/tailestim/estimators/pickands.py @@ -1,36 +1,40 @@ """Pickands estimator implementation for tail index estimation.""" + +from typing import Any, Dict, Tuple + import numpy as np -from typing import Dict, Any, Tuple -from .result import TailEstimatorResult + from .base import BaseTailEstimator +from .result import TailEstimatorResult from .tail_methods import pickands_estimator as pickands_estimate + class PickandsEstimator(BaseTailEstimator): """Pickands estimator for tail index estimation. - + This class implements the Pickands estimator, which is a simple method that does not use bootstrap procedures. Note that estimates can only be calculated up to the floor(n/4)-th order statistic. - + Parameters ---------- **kwargs : dict Additional parameters (not used by this estimator). """ - + def __init__(self, **kwargs): # Pickands estimator doesn't use bootstrap super().__init__(bootstrap=False, **kwargs) def _estimate(self, ordered_data: np.ndarray) -> Tuple: """Estimate the tail index using the Pickands estimator. - + Parameters ---------- ordered_data : np.ndarray Data array in decreasing order. - + Returns ------- Tuple @@ -40,19 +44,17 @@ def _estimate(self, ordered_data: np.ndarray) -> Tuple: def get_params(self) -> Dict[str, Any]: """Get the parameters of the estimator. - + Returns ------- dict Dictionary containing the parameters of the estimator. """ - return { - **self.kwargs - } + return {**self.kwargs} def get_result(self) -> TailEstimatorResult: """Get the estimated parameters. - + Attributes ---------- estimator : BaseTailEstimator @@ -68,13 +70,9 @@ def get_result(self) -> TailEstimatorResult: """ if self.results is None: raise ValueError("Model not fitted yet. Call fit() first.") - + k_arr, xi_arr = self.results - - res = { - 'estimator': self, - 'k_arr_': k_arr, - 'xi_arr_': xi_arr - } - - return TailEstimatorResult(res) \ No newline at end of file + + res = {"estimator": self, "k_arr_": k_arr, "xi_arr_": xi_arr} + + return TailEstimatorResult(res) diff --git a/src/tailestim/estimators/plot/plot_methods.py b/src/tailestim/estimators/plot/plot_methods.py index e86cb81..b646ad7 100644 --- a/src/tailestim/estimators/plot/plot_methods.py +++ b/src/tailestim/estimators/plot/plot_methods.py @@ -1,13 +1,25 @@ import os + import numpy as np from matplotlib import pyplot as plt -def make_plots(ordered_data, results, output_file_path, alpha=0.6, - bootstrap_flag=True, diagn_plots=False, theta1=0.01, theta2=0.99, - verbose=False, noise_flag=True, savedata=False): + +def make_plots( + ordered_data, + results, + output_file_path, + alpha=0.6, + bootstrap_flag=True, + diagn_plots=False, + theta1=0.01, + theta2=0.99, + verbose=False, + noise_flag=True, + savedata=False, +): """ Create plots for tail index estimation. - + Parameters ---------- ordered_data : numpy.ndarray @@ -17,7 +29,7 @@ def make_plots(ordered_data, results, output_file_path, alpha=0.6, output_file_path : str or None, optional File path to which plots should be saved. If None, the figure is not saved. alpha : float, optional - Parameter controlling the amount of "smoothing" for the kernel-type estimator. + Parameter controlling the amount of "smoothing" for the kernel-type estimator. Should be greater than 0.5, default is 0.6. bootstrap_flag : bool, optional Flag to switch on/off double-bootstrap procedure, default is True. @@ -36,7 +48,7 @@ def make_plots(ordered_data, results, output_file_path, alpha=0.6, data point. Used for integer-valued sequences with p = 1, default is True. savedata : bool, optional Flag to save data files in the directory with plots, default is False. - + Returns ------- fig : matplotlib.figure.Figure @@ -44,70 +56,97 @@ def make_plots(ordered_data, results, output_file_path, alpha=0.6, """ # Extract results for plotting - x_pdf, y_pdf = results['pdf']['x'], results['pdf']['y'] - x_ccdf, y_ccdf = results['ccdf']['x'], results['ccdf']['y'] - ordered_data = results['ordered_data'] + x_pdf, y_pdf = results["pdf"]["x"], results["pdf"]["y"] + x_ccdf, y_ccdf = results["ccdf"]["x"], results["ccdf"]["y"] + ordered_data = results["ordered_data"] if noise_flag: - discrete_ordered_data = results['discrete_ordered_data'] - + discrete_ordered_data = results["discrete_ordered_data"] + # Pickands results - k_p_arr, xi_p_arr = results['pickands']['k_arr_'], results['pickands']['xi_arr_'] - + k_p_arr, xi_p_arr = results["pickands"]["k_arr_"], results["pickands"]["xi_arr_"] + # Smooth Hill results - k_sh_arr, xi_sh_arr = results['smooth_hill']['k_arr_'], results['smooth_hill']['xi_arr_'] - + k_sh_arr, xi_sh_arr = ( + results["smooth_hill"]["k_arr_"], + results["smooth_hill"]["xi_arr_"], + ) + # Hill results - k_h_arr, xi_h_arr = results['hill']['k_arr_'], results['hill']['xi_arr_'] - if bootstrap_flag and 'k_star_' in results['hill']: - k_h_star = results['hill']['k_star_'] - xi_h_star = results['hill']['xi_star_'] - x1_h_arr = results['hill']['bootstrap_results_']['first_bootstrap_']['x_arr_'] - n1_h_amse = results['hill']['bootstrap_results_']['first_bootstrap_']['amse_'] - k1_h = results['hill']['bootstrap_results_']['first_bootstrap_']['k_min_'] - max_h_index1 = results['hill']['bootstrap_results_']['first_bootstrap_']['max_index_'] - x2_h_arr = results['hill']['bootstrap_results_']['second_bootstrap_']['x_arr_'] - n2_h_amse = results['hill']['bootstrap_results_']['second_bootstrap_']['amse_'] - k2_h = results['hill']['bootstrap_results_']['second_bootstrap_']['k_min_'] - max_h_index2 = results['hill']['bootstrap_results_']['second_bootstrap_']['max_index_'] + k_h_arr, xi_h_arr = results["hill"]["k_arr_"], results["hill"]["xi_arr_"] + if bootstrap_flag and "k_star_" in results["hill"]: + k_h_star = results["hill"]["k_star_"] + xi_h_star = results["hill"]["xi_star_"] + x1_h_arr = results["hill"]["bootstrap_results_"]["first_bootstrap_"]["x_arr_"] + n1_h_amse = results["hill"]["bootstrap_results_"]["first_bootstrap_"]["amse_"] + k1_h = results["hill"]["bootstrap_results_"]["first_bootstrap_"]["k_min_"] + max_h_index1 = results["hill"]["bootstrap_results_"]["first_bootstrap_"][ + "max_index_" + ] + x2_h_arr = results["hill"]["bootstrap_results_"]["second_bootstrap_"]["x_arr_"] + n2_h_amse = results["hill"]["bootstrap_results_"]["second_bootstrap_"]["amse_"] + k2_h = results["hill"]["bootstrap_results_"]["second_bootstrap_"]["k_min_"] + max_h_index2 = results["hill"]["bootstrap_results_"]["second_bootstrap_"][ + "max_index_" + ] else: k_h_star = None xi_h_star = None x1_h_arr, n1_h_amse, k1_h, max_h_index1 = None, None, None, None x2_h_arr, n2_h_amse, k2_h, max_h_index2 = None, None, None, None - + # Moments results - k_m_arr, xi_m_arr = results['moments']['k_arr_'], results['moments']['xi_arr_'] - if bootstrap_flag and 'k_star_' in results['moments']: - k_m_star = results['moments']['k_star_'] - xi_m_star = results['moments']['xi_star_'] - x1_m_arr = results['moments']['bootstrap_results_']['first_bootstrap_']['x_arr_'] - n1_m_amse = results['moments']['bootstrap_results_']['first_bootstrap_']['amse_'] - k1_m = results['moments']['bootstrap_results_']['first_bootstrap_']['k_min_'] - max_m_index1 = results['moments']['bootstrap_results_']['first_bootstrap_']['max_index_'] - x2_m_arr = results['moments']['bootstrap_results_']['second_bootstrap_']['x_arr_'] - n2_m_amse = results['moments']['bootstrap_results_']['second_bootstrap_']['amse_'] - k2_m = results['moments']['bootstrap_results_']['second_bootstrap_']['k_min_'] - max_m_index2 = results['moments']['bootstrap_results_']['second_bootstrap_']['max_index_'] + k_m_arr, xi_m_arr = results["moments"]["k_arr_"], results["moments"]["xi_arr_"] + if bootstrap_flag and "k_star_" in results["moments"]: + k_m_star = results["moments"]["k_star_"] + xi_m_star = results["moments"]["xi_star_"] + x1_m_arr = results["moments"]["bootstrap_results_"]["first_bootstrap_"][ + "x_arr_" + ] + n1_m_amse = results["moments"]["bootstrap_results_"]["first_bootstrap_"][ + "amse_" + ] + k1_m = results["moments"]["bootstrap_results_"]["first_bootstrap_"]["k_min_"] + max_m_index1 = results["moments"]["bootstrap_results_"]["first_bootstrap_"][ + "max_index_" + ] + x2_m_arr = results["moments"]["bootstrap_results_"]["second_bootstrap_"][ + "x_arr_" + ] + n2_m_amse = results["moments"]["bootstrap_results_"]["second_bootstrap_"][ + "amse_" + ] + k2_m = results["moments"]["bootstrap_results_"]["second_bootstrap_"]["k_min_"] + max_m_index2 = results["moments"]["bootstrap_results_"]["second_bootstrap_"][ + "max_index_" + ] else: k_m_star = None xi_m_star = None x1_m_arr, n1_m_amse, k1_m, max_m_index1 = None, None, None, None x2_m_arr, n2_m_amse, k2_m, max_m_index2 = None, None, None, None - + # Kernel results - k_k_arr, xi_k_arr = results['kernel']['k_arr_'], results['kernel']['xi_arr_'] - if bootstrap_flag and 'k_star_' in results['kernel']: - k_k_star = results['kernel']['k_star_'] - xi_k_star = results['kernel']['xi_star_'] - x1_k_arr = results['kernel']['bootstrap_results_']['first_bootstrap_']['x_arr_'] - n1_k_amse = results['kernel']['bootstrap_results_']['first_bootstrap_']['amse_'] - h1 = results['kernel']['bootstrap_results_']['first_bootstrap_']['h_min'] - max_k_index1 = results['kernel']['bootstrap_results_']['first_bootstrap_']['max_index_'] - x2_k_arr = results['kernel']['bootstrap_results_']['second_bootstrap_']['x_arr_'] - n2_k_amse = results['kernel']['bootstrap_results_']['second_bootstrap_']['amse_'] - h2 = results['kernel']['bootstrap_results_']['second_bootstrap_']['h_min'] - max_k_index2 = results['kernel']['bootstrap_results_']['second_bootstrap_']['max_index_'] - k_k1_star = results['kernel'].get('k_k1_star') + k_k_arr, xi_k_arr = results["kernel"]["k_arr_"], results["kernel"]["xi_arr_"] + if bootstrap_flag and "k_star_" in results["kernel"]: + k_k_star = results["kernel"]["k_star_"] + xi_k_star = results["kernel"]["xi_star_"] + x1_k_arr = results["kernel"]["bootstrap_results_"]["first_bootstrap_"]["x_arr_"] + n1_k_amse = results["kernel"]["bootstrap_results_"]["first_bootstrap_"]["amse_"] + h1 = results["kernel"]["bootstrap_results_"]["first_bootstrap_"]["h_min"] + max_k_index1 = results["kernel"]["bootstrap_results_"]["first_bootstrap_"][ + "max_index_" + ] + x2_k_arr = results["kernel"]["bootstrap_results_"]["second_bootstrap_"][ + "x_arr_" + ] + n2_k_amse = results["kernel"]["bootstrap_results_"]["second_bootstrap_"][ + "amse_" + ] + h2 = results["kernel"]["bootstrap_results_"]["second_bootstrap_"]["h_min"] + max_k_index2 = results["kernel"]["bootstrap_results_"]["second_bootstrap_"][ + "max_index_" + ] + k_k1_star = results["kernel"].get("k_k1_star") else: k_k_star = None xi_k_star = None @@ -118,77 +157,96 @@ def make_plots(ordered_data, results, output_file_path, alpha=0.6, if savedata: output_dir = os.path.dirname(os.path.realpath(output_file_path)) output_name = os.path.splitext(os.path.basename(output_file_path))[0] - + # Save PDF data - with open(os.path.join(output_dir+"/"+output_name+"_pdf.dat"), "w") as f: + with open(os.path.join(output_dir + "/" + output_name + "_pdf.dat"), "w") as f: for i in range(len(x_pdf)): f.write(str(x_pdf[i]) + " " + str(y_pdf[i]) + "\n") - + # Save CCDF data - with open(os.path.join(output_dir+"/"+output_name+"_ccdf.dat"), "w") as f: + with open(os.path.join(output_dir + "/" + output_name + "_ccdf.dat"), "w") as f: for i in range(len(x_ccdf)): f.write(str(x_ccdf[i]) + " " + str(y_ccdf[i]) + "\n") - + # Save additional data files if needed if savedata: # Save Pickands data - with open(os.path.join(output_dir+"/"+output_name+"_pickands.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_pickands.dat"), "w" + ) as f: for i in range(len(k_p_arr)): f.write(str(k_p_arr[i]) + " " + str(xi_p_arr[i]) + "\n") - + # Save smooth Hill data - with open(os.path.join(output_dir+"/"+output_name+"_sm_hill.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_sm_hill.dat"), "w" + ) as f: for i in range(len(k_sh_arr)): f.write(str(k_sh_arr[i]) + " " + str(xi_sh_arr[i]) + "\n") - + # Save Hill data - with open(os.path.join(output_dir+"/"+output_name+"_adj_hill_plot.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_adj_hill_plot.dat"), "w" + ) as f: for i in range(len(k_h_arr)): f.write(str(k_h_arr[i]) + " " + str(xi_h_arr[i]) + "\n") if bootstrap_flag and k_h_star is not None: - with open(os.path.join(output_dir+"/"+output_name+"_adj_hill_estimate.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_adj_hill_estimate.dat"), + "w", + ) as f: f.write(str(k_h_star) + " " + str(xi_h_star) + "\n") - + # Save moments data - with open(os.path.join(output_dir+"/"+output_name+"_mom_plot.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_mom_plot.dat"), "w" + ) as f: for i in range(len(k_m_arr)): f.write(str(k_m_arr[i]) + " " + str(xi_m_arr[i]) + "\n") if bootstrap_flag and k_m_star is not None: - with open(os.path.join(output_dir+"/"+output_name+"_mom_estimate.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_mom_estimate.dat"), "w" + ) as f: f.write(str(k_m_star) + " " + str(xi_m_star) + "\n") - + # Save kernel data - with open(os.path.join(output_dir+"/"+output_name+"_kern_plot.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_kern_plot.dat"), "w" + ) as f: for i in range(len(k_k_arr)): f.write(str(k_k_arr[i]) + " " + str(xi_k_arr[i]) + "\n") if bootstrap_flag and k_k1_star is not None: - with open(os.path.join(output_dir+"/"+output_name+"_kern_estimate.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_kern_estimate.dat"), "w" + ) as f: f.write(str(k_k_arr[k_k1_star]) + " " + str(xi_k_arr[k_k1_star]) + "\n") - + # plotting part if verbose: print("Making plots...") - fig, axes = plt.subplots(3, 2, figsize = (12, 16)) + fig, axes = plt.subplots(3, 2, figsize=(12, 16)) for ax in axes.reshape(-1): - ax.tick_params(direction='out', length=6, width=1.5, - labelsize = 12, which = 'major') - ax.tick_params(direction='out', length=3, width=1, which = 'minor') + ax.tick_params( + direction="out", length=6, width=1.5, labelsize=12, which="major" + ) + ax.tick_params(direction="out", length=3, width=1, which="minor") [i.set_linewidth(1.5) for i in ax.spines.values()] # plot PDF - axes[0,0].set_xlabel(r"Degree $k$", fontsize = 20) - axes[0,0].set_ylabel(r"$P(k)$", fontsize = 20) - axes[0,0].loglog(x_pdf, y_pdf, color = "#386cb0", marker = "s", - lw = 1.5, markeredgecolor = "black") + axes[0, 0].set_xlabel(r"Degree $k$", fontsize=20) + axes[0, 0].set_ylabel(r"$P(k)$", fontsize=20) + axes[0, 0].loglog( + x_pdf, y_pdf, color="#386cb0", marker="s", lw=1.5, markeredgecolor="black" + ) # plot CCDF - axes[0,1].set_xlabel(r"Degree $k$", fontsize = 20) - axes[0,1].set_ylabel(r"$\bar{F}(k)$", fontsize = 20) - axes[0,1].set_xscale("log") - axes[0,1].set_yscale("log") - axes[0,1].step(x_ccdf, y_ccdf, color = "#386cb0", lw = 1.5) - + axes[0, 1].set_xlabel(r"Degree $k$", fontsize=20) + axes[0, 1].set_ylabel(r"$\bar{F}(k)$", fontsize=20) + axes[0, 1].set_xscale("log") + axes[0, 1].set_yscale("log") + axes[0, 1].step(x_ccdf, y_ccdf, color="#386cb0", lw=1.5) + # draw scalings # Draw Hill scaling only if bootstrap is enabled and k_h_star is available if bootstrap_flag and k_h_star is not None: @@ -199,16 +257,29 @@ def make_plots(ordered_data, results, output_file_path, alpha=0.6, x = x_ccdf[np.where(x_ccdf >= xmin)] if len(x) > 0: # Make sure we have points to plot l = np.mean(y_ccdf[np.where(x == xmin)]) - alpha = 1./xi_h_star + alpha = 1.0 / xi_h_star if xi_h_star > 0: - axes[0,1].plot(x, [l*(float(xmin)/k)**alpha for k in x], - color = '#fb8072', ls = '--', lw = 2, - label = r"Adj. Hill Scaling $(\alpha="+\ - str(np.round(1./xi_h_star, decimals = 3))+r")$") - axes[0,1].plot((x[-1]), [l*(float(xmin)/x[-1])**(alpha)], - color = "#fb8072", ls = 'none', marker = 'o', - markerfacecolor = 'none', markeredgecolor = "#fb8072", - markeredgewidth = 3, markersize = 10) + axes[0, 1].plot( + x, + [l * (float(xmin) / k) ** alpha for k in x], + color="#fb8072", + ls="--", + lw=2, + label=r"Adj. Hill Scaling $(\alpha=" + + str(np.round(1.0 / xi_h_star, decimals=3)) + + r")$", + ) + axes[0, 1].plot( + (x[-1]), + [l * (float(xmin) / x[-1]) ** (alpha)], + color="#fb8072", + ls="none", + marker="o", + markerfacecolor="none", + markeredgecolor="#fb8072", + markeredgewidth=3, + markersize=10, + ) # Draw Moments scaling only if bootstrap is enabled and k_m_star is available if bootstrap_flag and k_m_star is not None: if noise_flag: @@ -218,214 +289,332 @@ def make_plots(ordered_data, results, output_file_path, alpha=0.6, x = x_ccdf[np.where(x_ccdf >= xmin)] if len(x) > 0: # Make sure we have points to plot l = np.mean(y_ccdf[np.where(x == xmin)]) - alpha = 1./xi_m_star + alpha = 1.0 / xi_m_star if xi_m_star > 0: - axes[0,1].plot(x, [l*(float(xmin)/k)**alpha for k in x], - color = '#8dd3c7', ls = '--', lw = 2, - label = r"Moments Scaling $(\alpha="+\ - str(np.round(1./xi_m_star, decimals = 3))+r")$") - axes[0,1].plot((x[-1]), [l*(float(xmin)/x[-1])**(alpha)], - color = "#8dd3c7", ls = 'none', marker = 'o', - markerfacecolor = 'none', markeredgecolor = "#8dd3c7", - markeredgewidth = 3, markersize = 10) + axes[0, 1].plot( + x, + [l * (float(xmin) / k) ** alpha for k in x], + color="#8dd3c7", + ls="--", + lw=2, + label=r"Moments Scaling $(\alpha=" + + str(np.round(1.0 / xi_m_star, decimals=3)) + + r")$", + ) + axes[0, 1].plot( + (x[-1]), + [l * (float(xmin) / x[-1]) ** (alpha)], + color="#8dd3c7", + ls="none", + marker="o", + markerfacecolor="none", + markeredgecolor="#8dd3c7", + markeredgewidth=3, + markersize=10, + ) # Draw Kernel scaling only if bootstrap is enabled and k_k_star is available if bootstrap_flag and k_k_star is not None: if noise_flag: xmin = discrete_ordered_data[k_k_star] else: xmin = ordered_data[k_k_star] - + x = x_ccdf[np.where(x_ccdf >= xmin)] if len(x) > 0: # Make sure we have points to plot l = np.mean(y_ccdf[np.where(x == xmin)]) - alpha = 1./xi_k_star + alpha = 1.0 / xi_k_star if xi_k_star > 0: - axes[0,1].plot(x, [l*(float(xmin)/k)**alpha for k in x], - color = '#fdb462', ls = '--', lw = 2, - label = r"Kernel Scaling $(\alpha="+\ - str(np.round(1./xi_k_star, decimals = 3))+r")$") - axes[0,1].plot((x[-1]), [l*(float(xmin)/x[-1])**(alpha)], - color = "#fdb462", ls = 'none', marker = 'o', - markerfacecolor = 'none', markeredgecolor = "#fdb462", - markeredgewidth = 3, markersize = 10) - axes[0,1].legend(loc = 'best') + axes[0, 1].plot( + x, + [l * (float(xmin) / k) ** alpha for k in x], + color="#fdb462", + ls="--", + lw=2, + label=r"Kernel Scaling $(\alpha=" + + str(np.round(1.0 / xi_k_star, decimals=3)) + + r")$", + ) + axes[0, 1].plot( + (x[-1]), + [l * (float(xmin) / x[-1]) ** (alpha)], + color="#fdb462", + ls="none", + marker="o", + markerfacecolor="none", + markeredgecolor="#fdb462", + markeredgewidth=3, + markersize=10, + ) + axes[0, 1].legend(loc="best") # define min and max order statistics to plot - min_k = int(np.ceil(len(k_h_arr)**theta1)) - 1 - max_k = int(np.floor(len(k_h_arr)**theta2)) - 1 + min_k = int(np.ceil(len(k_h_arr) ** theta1)) - 1 + max_k = int(np.floor(len(k_h_arr) ** theta2)) - 1 # check if estimators' values are not too off in these bounds min_k_index = (np.abs(k_sh_arr - min_k)).argmin() max_k_index = (np.abs(k_sh_arr - max_k)).argmin() - if (xi_sh_arr[min_k_index] <= -3 or xi_sh_arr[min_k_index] >= 3): + if xi_sh_arr[min_k_index] <= -3 or xi_sh_arr[min_k_index] >= 3: indices_to_plot_sh = np.where((xi_sh_arr <= 3) & (xi_sh_arr >= -3)) - elif (xi_sh_arr[max_k_index] <= -3 or xi_sh_arr[max_k_index] >= 3): + elif xi_sh_arr[max_k_index] <= -3 or xi_sh_arr[max_k_index] >= 3: indices_to_plot_sh = np.where((xi_sh_arr <= 3) & (xi_sh_arr >= -3)) else: indices_to_plot_sh = np.where((k_sh_arr <= max_k) & (k_sh_arr >= min_k)) - axes[1,0].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize = 20) - axes[1,0].set_ylabel(r"Estimated $\hat{\xi}$", fontsize = 20) + axes[1, 0].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize=20) + axes[1, 0].set_ylabel(r"Estimated $\hat{\xi}$", fontsize=20) # plot smooth Hill - - axes[1,0].plot(k_sh_arr[indices_to_plot_sh], xi_sh_arr[indices_to_plot_sh], - color = "#b3de69", alpha = 0.8, label = "Smooth Hill", - zorder = 10) + + axes[1, 0].plot( + k_sh_arr[indices_to_plot_sh], + xi_sh_arr[indices_to_plot_sh], + color="#b3de69", + alpha=0.8, + label="Smooth Hill", + zorder=10, + ) # plot adjusted Hill # check if estimators' values are not too off in these bounds - if (xi_h_arr[min_k-1] <= -3 or xi_h_arr[min_k-1] >= 3): + if xi_h_arr[min_k - 1] <= -3 or xi_h_arr[min_k - 1] >= 3: indices_to_plot_h = np.where((xi_h_arr <= 3) & (xi_h_arr >= -3)) - elif (xi_h_arr[max_k-1] <= -3 or xi_h_arr[max_k-1] >= 3): + elif xi_h_arr[max_k - 1] <= -3 or xi_h_arr[max_k - 1] >= 3: indices_to_plot_h = np.where((xi_h_arr <= 3) & (xi_h_arr >= -3)) else: indices_to_plot_h = np.where((k_h_arr <= max_k) & (k_h_arr >= min_k)) - axes[1,0].plot(k_h_arr[indices_to_plot_h], xi_h_arr[indices_to_plot_h], - color = "#fb8072", alpha = 0.8, label = "Adjusted Hill", - zorder = 10) + axes[1, 0].plot( + k_h_arr[indices_to_plot_h], + xi_h_arr[indices_to_plot_h], + color="#fb8072", + alpha=0.8, + label="Adjusted Hill", + zorder=10, + ) if bootstrap_flag and k_h_star is not None: # Find the index of k_h_star in k_h_arr for plotting k_h_idx = np.argmin(np.abs(k_h_arr - k_h_star)) - axes[1,0].scatter([k_h_arr[k_h_idx]], [xi_h_arr[k_h_idx]], - color = "#fb8072", marker = "*", s = 100, - edgecolor = "black", zorder = 20, - label = r"$\widehat{\xi}^{Hill}="\ - +str(np.round([xi_h_arr[k_h_idx]][0], decimals = 3))\ - +r"$") - axes[1,0].legend(loc = "best") - - - axes[1,1].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize = 20) - axes[1,1].set_ylabel(r"Estimated $\hat{\xi}$", fontsize = 20) - axes[1,1].set_xscale("log") - + axes[1, 0].scatter( + [k_h_arr[k_h_idx]], + [xi_h_arr[k_h_idx]], + color="#fb8072", + marker="*", + s=100, + edgecolor="black", + zorder=20, + label=r"$\widehat{\xi}^{Hill}=" + + str(np.round([xi_h_arr[k_h_idx]][0], decimals=3)) + + r"$", + ) + axes[1, 0].legend(loc="best") + + axes[1, 1].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize=20) + axes[1, 1].set_ylabel(r"Estimated $\hat{\xi}$", fontsize=20) + axes[1, 1].set_xscale("log") + # plot smooth Hill - axes[1,1].plot(k_sh_arr[indices_to_plot_sh], xi_sh_arr[indices_to_plot_sh], - color = "#b3de69", alpha = 0.8, label = "Smooth Hill", - zorder = 10) + axes[1, 1].plot( + k_sh_arr[indices_to_plot_sh], + xi_sh_arr[indices_to_plot_sh], + color="#b3de69", + alpha=0.8, + label="Smooth Hill", + zorder=10, + ) # plot adjusted Hill indices_to_plot = np.where((k_h_arr <= max_k) & (k_h_arr >= min_k)) - axes[1,1].plot(k_h_arr[indices_to_plot_h], xi_h_arr[indices_to_plot_h], - color = "#fb8072", alpha = 0.8, label = "Adjusted Hill", - zorder = 10) + axes[1, 1].plot( + k_h_arr[indices_to_plot_h], + xi_h_arr[indices_to_plot_h], + color="#fb8072", + alpha=0.8, + label="Adjusted Hill", + zorder=10, + ) if bootstrap_flag and k_h_star is not None: # Find the index of k_h_star in k_h_arr for plotting k_h_idx = np.argmin(np.abs(k_h_arr - k_h_star)) - axes[1,1].scatter([k_h_arr[k_h_idx]], [xi_h_arr[k_h_idx]], - color = "#fb8072", marker = "*", s = 100, - edgecolor = "black", zorder = 20, - label = r"$\widehat{\xi}^{Hill}="\ - +str(np.round([xi_h_arr[k_h_idx]][0], decimals = 3))\ - +r"$") - axes[1,1].legend(loc = "best") - - axes[2,0].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize = 20) - axes[2,0].set_ylabel(r"Estimated $\hat{\xi}$", fontsize = 20) - #plot Pickands + axes[1, 1].scatter( + [k_h_arr[k_h_idx]], + [xi_h_arr[k_h_idx]], + color="#fb8072", + marker="*", + s=100, + edgecolor="black", + zorder=20, + label=r"$\widehat{\xi}^{Hill}=" + + str(np.round([xi_h_arr[k_h_idx]][0], decimals=3)) + + r"$", + ) + axes[1, 1].legend(loc="best") + + axes[2, 0].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize=20) + axes[2, 0].set_ylabel(r"Estimated $\hat{\xi}$", fontsize=20) + # plot Pickands min_k_index = (np.abs(k_p_arr - min_k)).argmin() max_k_index = (np.abs(k_p_arr - max_k)).argmin() - if (xi_p_arr[min_k_index] <= -3 or xi_p_arr[min_k_index] >= 3): + if xi_p_arr[min_k_index] <= -3 or xi_p_arr[min_k_index] >= 3: indices_to_plot_p = np.where((xi_p_arr <= 3) & (xi_p_arr >= -3)) - elif (xi_p_arr[max_k_index] <= -3 or xi_p_arr[max_k_index] >= 3): + elif xi_p_arr[max_k_index] <= -3 or xi_p_arr[max_k_index] >= 3: indices_to_plot_p = np.where((xi_p_arr <= 3) & (xi_p_arr >= -3)) else: indices_to_plot_p = np.where((k_p_arr <= max_k) & (k_p_arr >= min_k)) - axes[2,0].plot(k_p_arr[indices_to_plot_p], xi_p_arr[indices_to_plot_p], - color = "#bc80bd", alpha = 0.8, label = "Pickands", - zorder = 10) - #plot moments - if (xi_m_arr[min_k-1] <= -3 or xi_m_arr[min_k-1] >= 3): + axes[2, 0].plot( + k_p_arr[indices_to_plot_p], + xi_p_arr[indices_to_plot_p], + color="#bc80bd", + alpha=0.8, + label="Pickands", + zorder=10, + ) + # plot moments + if xi_m_arr[min_k - 1] <= -3 or xi_m_arr[min_k - 1] >= 3: indices_to_plot_m = np.where((xi_m_arr <= 3) & (xi_m_arr >= -3)) - elif (xi_m_arr[max_k-1] <= -3 or xi_m_arr[max_k-1] >= 3): + elif xi_m_arr[max_k - 1] <= -3 or xi_m_arr[max_k - 1] >= 3: indices_to_plot_m = np.where((xi_m_arr <= 3) & (xi_m_arr >= -3)) else: indices_to_plot_m = np.where((k_m_arr <= max_k) & (k_m_arr >= min_k)) - - axes[2,0].plot(k_m_arr[indices_to_plot_m], xi_m_arr[indices_to_plot_m], - color = "#8dd3c7", alpha = 0.8, label = "Moments", - zorder = 10) + + axes[2, 0].plot( + k_m_arr[indices_to_plot_m], + xi_m_arr[indices_to_plot_m], + color="#8dd3c7", + alpha=0.8, + label="Moments", + zorder=10, + ) if bootstrap_flag and k_m_star is not None: # Find the index of k_m_star in k_m_arr for plotting k_m_idx = np.argmin(np.abs(k_m_arr - k_m_star)) - axes[2,0].scatter([k_m_arr[k_m_idx]], [xi_m_arr[k_m_idx]], - color = "#8dd3c7", marker = "*", s = 100, - edgecolor = "black", zorder = 20, - label = r"$\widehat{\xi}^{Moments}="\ - +str(np.round([xi_m_arr[k_m_idx]][0], decimals = 3))\ - +r"$") - #plot kernel-type + axes[2, 0].scatter( + [k_m_arr[k_m_idx]], + [xi_m_arr[k_m_idx]], + color="#8dd3c7", + marker="*", + s=100, + edgecolor="black", + zorder=20, + label=r"$\widehat{\xi}^{Moments}=" + + str(np.round([xi_m_arr[k_m_idx]][0], decimals=3)) + + r"$", + ) + # plot kernel-type min_k_index = (np.abs(k_k_arr - min_k)).argmin() max_k_index = (np.abs(k_k_arr - max_k)).argmin() - if (xi_k_arr[min_k_index] <= -3 or xi_k_arr[min_k_index] >= 3): + if xi_k_arr[min_k_index] <= -3 or xi_k_arr[min_k_index] >= 3: indices_to_plot_k = np.where((xi_k_arr <= 3) & (xi_k_arr >= -3)) - elif (xi_k_arr[max_k_index] <= -3 or xi_k_arr[max_k_index] >= 3): + elif xi_k_arr[max_k_index] <= -3 or xi_k_arr[max_k_index] >= 3: indices_to_plot_k = np.where((xi_k_arr <= 3) & (xi_k_arr >= -3)) else: indices_to_plot_k = list(range(min_k_index, max_k_index)) - #indices_to_plot_k = np.where((xi_k_arr <= 3) & (xi_k_arr >= -3)) - axes[2,0].plot(k_k_arr[indices_to_plot_k], xi_k_arr[indices_to_plot_k], - color = "#fdb462", alpha = 0.8, label = "Kernel", - zorder = 10) + # indices_to_plot_k = np.where((xi_k_arr <= 3) & (xi_k_arr >= -3)) + axes[2, 0].plot( + k_k_arr[indices_to_plot_k], + xi_k_arr[indices_to_plot_k], + color="#fdb462", + alpha=0.8, + label="Kernel", + zorder=10, + ) if bootstrap_flag and k_k1_star is not None: - axes[2,0].scatter([k_k_arr[k_k1_star]], [xi_k_arr[k_k1_star]], - color = "#fdb462", marker = "*", s = 100, - edgecolor = "black", zorder = 20, - label = r"$\widehat{\xi}^{Kernel}="\ - +str(np.round([xi_k_arr[k_k1_star]][0], decimals = 3))\ - +r"$") - axes[2,0].legend(loc = "best") + axes[2, 0].scatter( + [k_k_arr[k_k1_star]], + [xi_k_arr[k_k1_star]], + color="#fdb462", + marker="*", + s=100, + edgecolor="black", + zorder=20, + label=r"$\widehat{\xi}^{Kernel}=" + + str(np.round([xi_k_arr[k_k1_star]][0], decimals=3)) + + r"$", + ) + axes[2, 0].legend(loc="best") # for clarity purposes, display only xi region between -1 and 1 - axes[2,0].set_ylim((-0.5,1.5)) - - axes[2,1].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize = 20) - axes[2,1].set_ylabel(r"Estimated $\hat{\xi}$", fontsize = 20) - axes[2,1].set_xscale("log") - - #plot Pickands - axes[2,1].plot(k_p_arr[indices_to_plot_p], xi_p_arr[indices_to_plot_p], - color = "#bc80bd", alpha = 0.8, label = "Pickands", - zorder = 10) - #plot moments - axes[2,1].plot(k_m_arr[indices_to_plot_m], xi_m_arr[indices_to_plot_m], - color = "#8dd3c7", alpha = 0.8, label = "Moments", - zorder = 10) + axes[2, 0].set_ylim((-0.5, 1.5)) + + axes[2, 1].set_xlabel(r"Number of Order Statistics $\kappa$", fontsize=20) + axes[2, 1].set_ylabel(r"Estimated $\hat{\xi}$", fontsize=20) + axes[2, 1].set_xscale("log") + + # plot Pickands + axes[2, 1].plot( + k_p_arr[indices_to_plot_p], + xi_p_arr[indices_to_plot_p], + color="#bc80bd", + alpha=0.8, + label="Pickands", + zorder=10, + ) + # plot moments + axes[2, 1].plot( + k_m_arr[indices_to_plot_m], + xi_m_arr[indices_to_plot_m], + color="#8dd3c7", + alpha=0.8, + label="Moments", + zorder=10, + ) if bootstrap_flag and k_m_star is not None: # Find the index of k_m_star in k_m_arr for plotting k_m_idx = np.argmin(np.abs(k_m_arr - k_m_star)) - axes[2,1].scatter([k_m_arr[k_m_idx]], [xi_m_arr[k_m_idx]], - color = "#8dd3c7", marker = "*", s = 100, - edgecolor = "black", zorder = 20, - label = r"$\widehat{\xi}^{Moments}="\ - +str(np.round([xi_m_arr[k_m_idx]][0], decimals = 3))\ - +r"$") - #plot kernel-type - axes[2,1].plot(k_k_arr[indices_to_plot_k], xi_k_arr[indices_to_plot_k], - color = "#fdb462", alpha = 0.8, label = "Kernel", - zorder = 10) + axes[2, 1].scatter( + [k_m_arr[k_m_idx]], + [xi_m_arr[k_m_idx]], + color="#8dd3c7", + marker="*", + s=100, + edgecolor="black", + zorder=20, + label=r"$\widehat{\xi}^{Moments}=" + + str(np.round([xi_m_arr[k_m_idx]][0], decimals=3)) + + r"$", + ) + # plot kernel-type + axes[2, 1].plot( + k_k_arr[indices_to_plot_k], + xi_k_arr[indices_to_plot_k], + color="#fdb462", + alpha=0.8, + label="Kernel", + zorder=10, + ) if bootstrap_flag and k_k1_star is not None: - axes[2,1].scatter([k_k_arr[k_k1_star]], [xi_k_arr[k_k1_star]], - color = "#fdb462", marker = "*", s = 100, - edgecolor = "black", zorder = 20, - label = r"$\widehat{\xi}^{Kernel}="\ - +str(np.round([xi_k_arr[k_k1_star]][0], decimals = 3))\ - +r"$") + axes[2, 1].scatter( + [k_k_arr[k_k1_star]], + [xi_k_arr[k_k1_star]], + color="#fdb462", + marker="*", + s=100, + edgecolor="black", + zorder=20, + label=r"$\widehat{\xi}^{Kernel}=" + + str(np.round([xi_k_arr[k_k1_star]][0], decimals=3)) + + r"$", + ) # for clarity purposes, display only xi region between -1 and 1 - axes[2,1].set_ylim((-0.5,1.5)) - axes[2,1].legend(loc = "best") + axes[2, 1].set_ylim((-0.5, 1.5)) + axes[2, 1].legend(loc="best") # If diagnostic plots are requested, they are now handled by make_diagnostic_plots() if diagn_plots and bootstrap_flag and verbose: print("Diagnostic plots will be created separately using plot_diagnostics()") - - fig.tight_layout(pad = 0.2) + + fig.tight_layout(pad=0.2) if output_file_path is not None: fig.savefig(output_file_path) - + return fig, axes -def make_diagnostic_plots(results, output_file_path=None, - hsteps=200, bootstrap_flag=True, verbose=False, noise_flag=True, - savedata=False): + +def make_diagnostic_plots( + results, + output_file_path=None, + hsteps=200, + bootstrap_flag=True, + verbose=False, + noise_flag=True, + savedata=False, +): """ Create diagnostic plots for tail index estimation. - + Parameters ---------- results : dict @@ -443,7 +632,7 @@ def make_diagnostic_plots(results, output_file_path=None, data point. Used for integer-valued sequences with p = 1, default is True. savedata : bool, optional Flag to save data files in the directory with plots, default is False. - + Returns ------- fig : matplotlib.figure.Figure or None @@ -454,81 +643,105 @@ def make_diagnostic_plots(results, output_file_path=None, if verbose: print("Diagnostic plots require bootstrap to be enabled.") return None, None - + # Extract results for plotting - ordered_data = results['ordered_data'] - if noise_flag and 'discrete_ordered_data' in results: - discrete_ordered_data = results['discrete_ordered_data'] - + ordered_data = results["ordered_data"] + if noise_flag and "discrete_ordered_data" in results: + discrete_ordered_data = results["discrete_ordered_data"] + # Hill results - k_h_arr, xi_h_arr = results['hill']['k_arr_'], results['hill']['xi_arr_'] - if 'k_star_' in results['hill']: - k_h_star = results['hill']['k_star_'] - xi_h_star = results['hill']['xi_star_'] - x1_h_arr = results['hill']['bootstrap_results_']['first_bootstrap_']['x_arr_'] - n1_h_amse = results['hill']['bootstrap_results_']['first_bootstrap_']['amse_'] - k1_h = results['hill']['bootstrap_results_']['first_bootstrap_']['k_min_'] - max_h_index1 = results['hill']['bootstrap_results_']['first_bootstrap_']['max_index_'] - x2_h_arr = results['hill']['bootstrap_results_']['second_bootstrap_']['x_arr_'] - n2_h_amse = results['hill']['bootstrap_results_']['second_bootstrap_']['amse_'] - k2_h = results['hill']['bootstrap_results_']['second_bootstrap_']['k_min_'] - max_h_index2 = results['hill']['bootstrap_results_']['second_bootstrap_']['max_index_'] + k_h_arr, xi_h_arr = results["hill"]["k_arr_"], results["hill"]["xi_arr_"] + if "k_star_" in results["hill"]: + k_h_star = results["hill"]["k_star_"] + xi_h_star = results["hill"]["xi_star_"] + x1_h_arr = results["hill"]["bootstrap_results_"]["first_bootstrap_"]["x_arr_"] + n1_h_amse = results["hill"]["bootstrap_results_"]["first_bootstrap_"]["amse_"] + k1_h = results["hill"]["bootstrap_results_"]["first_bootstrap_"]["k_min_"] + max_h_index1 = results["hill"]["bootstrap_results_"]["first_bootstrap_"][ + "max_index_" + ] + x2_h_arr = results["hill"]["bootstrap_results_"]["second_bootstrap_"]["x_arr_"] + n2_h_amse = results["hill"]["bootstrap_results_"]["second_bootstrap_"]["amse_"] + k2_h = results["hill"]["bootstrap_results_"]["second_bootstrap_"]["k_min_"] + max_h_index2 = results["hill"]["bootstrap_results_"]["second_bootstrap_"][ + "max_index_" + ] else: k_h_star = None xi_h_star = None x1_h_arr, n1_h_amse, k1_h, max_h_index1 = None, None, None, None x2_h_arr, n2_h_amse, k2_h, max_h_index2 = None, None, None, None - + # Moments results - k_m_arr, xi_m_arr = results['moments']['k_arr_'], results['moments']['xi_arr_'] - if 'k_star_' in results['moments']: - k_m_star = results['moments']['k_star_'] - xi_m_star = results['moments']['xi_star_'] - x1_m_arr = results['moments']['bootstrap_results_']['first_bootstrap_']['x_arr_'] - n1_m_amse = results['moments']['bootstrap_results_']['first_bootstrap_']['amse_'] - k1_m = results['moments']['bootstrap_results_']['first_bootstrap_']['k_min_'] - max_m_index1 = results['moments']['bootstrap_results_']['first_bootstrap_']['max_index_'] - x2_m_arr = results['moments']['bootstrap_results_']['second_bootstrap_']['x_arr_'] - n2_m_amse = results['moments']['bootstrap_results_']['second_bootstrap_']['amse_'] - k2_m = results['moments']['bootstrap_results_']['second_bootstrap_']['k_min_'] - max_m_index2 = results['moments']['bootstrap_results_']['second_bootstrap_']['max_index_'] + k_m_arr, xi_m_arr = results["moments"]["k_arr_"], results["moments"]["xi_arr_"] + if "k_star_" in results["moments"]: + k_m_star = results["moments"]["k_star_"] + xi_m_star = results["moments"]["xi_star_"] + x1_m_arr = results["moments"]["bootstrap_results_"]["first_bootstrap_"][ + "x_arr_" + ] + n1_m_amse = results["moments"]["bootstrap_results_"]["first_bootstrap_"][ + "amse_" + ] + k1_m = results["moments"]["bootstrap_results_"]["first_bootstrap_"]["k_min_"] + max_m_index1 = results["moments"]["bootstrap_results_"]["first_bootstrap_"][ + "max_index_" + ] + x2_m_arr = results["moments"]["bootstrap_results_"]["second_bootstrap_"][ + "x_arr_" + ] + n2_m_amse = results["moments"]["bootstrap_results_"]["second_bootstrap_"][ + "amse_" + ] + k2_m = results["moments"]["bootstrap_results_"]["second_bootstrap_"]["k_min_"] + max_m_index2 = results["moments"]["bootstrap_results_"]["second_bootstrap_"][ + "max_index_" + ] else: k_m_star = None xi_m_star = None x1_m_arr, n1_m_amse, k1_m, max_m_index1 = None, None, None, None x2_m_arr, n2_m_amse, k2_m, max_m_index2 = None, None, None, None - + # Kernel results - k_k_arr, xi_k_arr = results['kernel']['k_arr_'], results['kernel']['xi_arr_'] - if 'k_star_' in results['kernel']: - k_k_star = results['kernel']['k_star_'] - xi_k_star = results['kernel']['xi_star_'] - x1_k_arr = results['kernel']['bootstrap_results_']['first_bootstrap_']['x_arr_'] - n1_k_amse = results['kernel']['bootstrap_results_']['first_bootstrap_']['amse_'] - h1 = results['kernel']['bootstrap_results_']['first_bootstrap_']['h_min'] - max_k_index1 = results['kernel']['bootstrap_results_']['first_bootstrap_']['max_index_'] - x2_k_arr = results['kernel']['bootstrap_results_']['second_bootstrap_']['x_arr_'] - n2_k_amse = results['kernel']['bootstrap_results_']['second_bootstrap_']['amse_'] - h2 = results['kernel']['bootstrap_results_']['second_bootstrap_']['h_min'] - max_k_index2 = results['kernel']['bootstrap_results_']['second_bootstrap_']['max_index_'] - k_k1_star = results['kernel'].get('k_k1_star') + k_k_arr, xi_k_arr = results["kernel"]["k_arr_"], results["kernel"]["xi_arr_"] + if "k_star_" in results["kernel"]: + k_k_star = results["kernel"]["k_star_"] + xi_k_star = results["kernel"]["xi_star_"] + x1_k_arr = results["kernel"]["bootstrap_results_"]["first_bootstrap_"]["x_arr_"] + n1_k_amse = results["kernel"]["bootstrap_results_"]["first_bootstrap_"]["amse_"] + h1 = results["kernel"]["bootstrap_results_"]["first_bootstrap_"]["h_min"] + max_k_index1 = results["kernel"]["bootstrap_results_"]["first_bootstrap_"][ + "max_index_" + ] + x2_k_arr = results["kernel"]["bootstrap_results_"]["second_bootstrap_"][ + "x_arr_" + ] + n2_k_amse = results["kernel"]["bootstrap_results_"]["second_bootstrap_"][ + "amse_" + ] + h2 = results["kernel"]["bootstrap_results_"]["second_bootstrap_"]["h_min"] + max_k_index2 = results["kernel"]["bootstrap_results_"]["second_bootstrap_"][ + "max_index_" + ] + k_k1_star = results["kernel"].get("k_k1_star") else: k_k_star = None xi_k_star = None x1_k_arr, n1_k_amse, h1, max_k_index1 = None, None, None, None x2_k_arr, n2_k_amse, h2, max_k_index2 = None, None, None, None k_k1_star = None - + # Setup for saving data if needed if savedata and output_file_path is not None: output_dir = os.path.dirname(os.path.realpath(output_file_path)) output_name = os.path.splitext(os.path.basename(output_file_path))[0] - + # Create diagnostic plots if verbose: print("Creating diagnostic plots...") - - fig_d, axes_d = plt.subplots(1, 3, figsize = (18, 6)) + + fig_d, axes_d = plt.subplots(1, 3, figsize=(18, 6)) # filter out boundary values using theta parameters for Hill if x1_h_arr is not None and x2_h_arr is not None: @@ -536,163 +749,356 @@ def make_diagnostic_plots(results, output_file_path=None, max_k1 = len(x1_h_arr) - 1 min_k2 = 2 max_k2 = len(x2_h_arr) - 1 - + axes_d[0].set_yscale("log") axes_d[0].set_xscale("log") axes_d[1].set_xscale("log") axes_d[2].set_xscale("log") - + if n1_h_amse is not None and x1_h_arr is not None: - n1_h_amse[np.where((n1_h_amse == np.inf) |\ - (n1_h_amse == -np.inf))] = np.nan - axes_d[0].set_ylim((0.1*np.nanmin(n1_h_amse[min_k1:max_k1]), 1.0)) - - axes_d[0].set_xlabel("Fraction of Bootstrap Order Statistics", - fontsize = 20) - axes_d[0].set_ylabel(r"$\langle AMSE \rangle$", fontsize = 20) - axes_d[0].set_title("Adjusted Hill Estimator", fontsize = 20) - + n1_h_amse[np.where((n1_h_amse == np.inf) | (n1_h_amse == -np.inf))] = np.nan + axes_d[0].set_ylim((0.1 * np.nanmin(n1_h_amse[min_k1:max_k1]), 1.0)) + + axes_d[0].set_xlabel("Fraction of Bootstrap Order Statistics", fontsize=20) + axes_d[0].set_ylabel(r"$\langle AMSE \rangle$", fontsize=20) + axes_d[0].set_title("Adjusted Hill Estimator", fontsize=20) + # plot AMSE and corresponding minimum for Hill if x1_h_arr is not None and n1_h_amse is not None and k1_h is not None: - axes_d[0].plot(x1_h_arr[min_k1:max_k1], n1_h_amse[min_k1:max_k1], - alpha = 0.5, lw = 1.5, - color = "#d55e00", label = r"$n_1$ samples") - axes_d[0].scatter([k1_h], [n1_h_amse[int(len(x1_h_arr)*k1_h)-1]], - color = "#d55e00", - marker = 'o', edgecolor = "black", alpha = 0.5, - label = r"Min for $n_1$ sample") - + axes_d[0].plot( + x1_h_arr[min_k1:max_k1], + n1_h_amse[min_k1:max_k1], + alpha=0.5, + lw=1.5, + color="#d55e00", + label=r"$n_1$ samples", + ) + axes_d[0].scatter( + [k1_h], + [n1_h_amse[int(len(x1_h_arr) * k1_h) - 1]], + color="#d55e00", + marker="o", + edgecolor="black", + alpha=0.5, + label=r"Min for $n_1$ sample", + ) + if x2_h_arr is not None and n2_h_amse is not None and k2_h is not None: - axes_d[0].plot(x2_h_arr[min_k2:max_k2], n2_h_amse[min_k2:max_k2], - alpha = 0.5, lw = 1.5, - color = "#0072b2", label = r"$n_2$ samples") - axes_d[0].scatter([k2_h], [n2_h_amse[int(len(x2_h_arr)*k2_h)-1]], - color = "#0072b2", - marker = 'o', edgecolor = "black", alpha = 0.5, - label = r"Min for $n_2$ sample") - + axes_d[0].plot( + x2_h_arr[min_k2:max_k2], + n2_h_amse[min_k2:max_k2], + alpha=0.5, + lw=1.5, + color="#0072b2", + label=r"$n_2$ samples", + ) + axes_d[0].scatter( + [k2_h], + [n2_h_amse[int(len(x2_h_arr) * k2_h) - 1]], + color="#0072b2", + marker="o", + edgecolor="black", + alpha=0.5, + label=r"Min for $n_2$ sample", + ) + if max_h_index1 is not None and x1_h_arr is not None: - axes_d[0].axvline(max_h_index1/float(len(x1_h_arr)), color = "#d55e00", - ls = '--', alpha = 0.5, - label = r"Minimization boundary for $n_1$ sample") - + axes_d[0].axvline( + max_h_index1 / float(len(x1_h_arr)), + color="#d55e00", + ls="--", + alpha=0.5, + label=r"Minimization boundary for $n_1$ sample", + ) + if max_h_index2 is not None and x2_h_arr is not None: - axes_d[0].axvline(max_h_index2/float(len(x2_h_arr)), color = "#0072b2", - ls = '--', alpha = 0.5, - label = r"Minimization boundary for $n_2$ sample") - - axes_d[0].legend(loc = "best") + axes_d[0].axvline( + max_h_index2 / float(len(x2_h_arr)), + color="#0072b2", + ls="--", + alpha=0.5, + label=r"Minimization boundary for $n_2$ sample", + ) + + axes_d[0].legend(loc="best") if savedata and output_file_path is not None: - with open(os.path.join(output_dir+"/"+output_name+"_adjhill_diagn1.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_adjhill_diagn1.dat"), "w" + ) as f: for i in range(len(x1_h_arr[min_k1:max_k1])): - f.write(str(x1_h_arr[min_k1:max_k1][i]) + " " + str(n1_h_amse[min_k1:max_k1][i]) + "\n") - with open(os.path.join(output_dir+"/"+output_name+"_adjhill_diagn2.dat"), "w") as f: + f.write( + str(x1_h_arr[min_k1:max_k1][i]) + + " " + + str(n1_h_amse[min_k1:max_k1][i]) + + "\n" + ) + with open( + os.path.join(output_dir + "/" + output_name + "_adjhill_diagn2.dat"), "w" + ) as f: for i in range(len(x2_h_arr[min_k2:max_k2])): - f.write(str(x2_h_arr[min_k2:max_k2][i]) + " " + str(n2_h_amse[min_k2:max_k2][i]) + "\n") - with open(os.path.join(output_dir+"/"+output_name+"_adjhill_diagn_points.dat"), "w") as f: - f.write("Min for n1 sample: "+str(k1_h)+" "+str(n1_h_amse[int(len(x1_h_arr)*k1_h)-1])+"\n") - f.write("Min for n2 sample: "+str(k2_h)+" "+str(n2_h_amse[int(len(x2_h_arr)*k2_h)-1])+"\n") - f.write("Minimization boundary for n1 sample: "+str(max_h_index1/float(len(x1_h_arr)))+"\n") - f.write("Minimization boundary for n2 sample: "+str(max_h_index2/float(len(x2_h_arr)))+"\n") + f.write( + str(x2_h_arr[min_k2:max_k2][i]) + + " " + + str(n2_h_amse[min_k2:max_k2][i]) + + "\n" + ) + with open( + os.path.join(output_dir + "/" + output_name + "_adjhill_diagn_points.dat"), + "w", + ) as f: + f.write( + "Min for n1 sample: " + + str(k1_h) + + " " + + str(n1_h_amse[int(len(x1_h_arr) * k1_h) - 1]) + + "\n" + ) + f.write( + "Min for n2 sample: " + + str(k2_h) + + " " + + str(n2_h_amse[int(len(x2_h_arr) * k2_h) - 1]) + + "\n" + ) + f.write( + "Minimization boundary for n1 sample: " + + str(max_h_index1 / float(len(x1_h_arr))) + + "\n" + ) + f.write( + "Minimization boundary for n2 sample: " + + str(max_h_index2 / float(len(x2_h_arr))) + + "\n" + ) # filter out boundary values using theta parameters for moments min_k1 = 2 max_k1 = len(x1_m_arr) - 1 min_k2 = 2 max_k2 = len(x2_m_arr) - 1 - n1_m_amse[np.where((n1_m_amse == np.inf) |\ - (n1_m_amse == -np.inf))] = np.nan + n1_m_amse[np.where((n1_m_amse == np.inf) | (n1_m_amse == -np.inf))] = np.nan axes_d[1].set_yscale("log") - axes_d[1].set_ylim((0.1*np.nanmin(n1_m_amse[min_k1:max_k1]), 1.0)) - axes_d[1].set_xlabel("Fraction of Bootstrap Order Statistics", - fontsize = 20) - axes_d[1].set_ylabel(r"$\langle AMSE \rangle$", fontsize = 20) - axes_d[1].set_title("Moments Estimator", fontsize = 20) + axes_d[1].set_ylim((0.1 * np.nanmin(n1_m_amse[min_k1:max_k1]), 1.0)) + axes_d[1].set_xlabel("Fraction of Bootstrap Order Statistics", fontsize=20) + axes_d[1].set_ylabel(r"$\langle AMSE \rangle$", fontsize=20) + axes_d[1].set_title("Moments Estimator", fontsize=20) # plot AMSE and corresponding minimum - axes_d[1].plot(x1_m_arr[min_k1:max_k1], n1_m_amse[min_k1:max_k1], - alpha = 0.5, lw = 1.5, - color = "#d55e00", label = r"$n_1$ samples") - axes_d[1].scatter([k1_m], [n1_m_amse[int(len(x1_m_arr)*k1_m)-1]], - color = "#d55e00", - marker = 'o', edgecolor = "black", alpha = 0.5, - label = r"Min for $n_1$ sample") - axes_d[1].plot(x2_m_arr[min_k2:max_k2], n2_m_amse[min_k2:max_k2], - alpha = 0.5, lw = 1.5, - color = "#0072b2", label = r"$n_2$ samples") - axes_d[1].scatter([k2_m], [n2_m_amse[int(len(x2_m_arr)*k2_m)-1]], - color = "#0072b2", - marker = 'o', edgecolor = "black", alpha = 0.5, - label = r"Min for $n_2$ sample") - axes_d[1].axvline(max_m_index1/float(len(x1_m_arr)), color = "#d55e00", - ls = '--', alpha = 0.5, - label = r"Minimization boundary for $n_1$ sample") - axes_d[1].axvline(max_m_index2/float(len(x2_m_arr)), color = "#0072b2", - ls = '--', alpha = 0.5, - label = r"Minimization boundary for $n_2$ sample") - axes_d[1].legend(loc = "best") + axes_d[1].plot( + x1_m_arr[min_k1:max_k1], + n1_m_amse[min_k1:max_k1], + alpha=0.5, + lw=1.5, + color="#d55e00", + label=r"$n_1$ samples", + ) + axes_d[1].scatter( + [k1_m], + [n1_m_amse[int(len(x1_m_arr) * k1_m) - 1]], + color="#d55e00", + marker="o", + edgecolor="black", + alpha=0.5, + label=r"Min for $n_1$ sample", + ) + axes_d[1].plot( + x2_m_arr[min_k2:max_k2], + n2_m_amse[min_k2:max_k2], + alpha=0.5, + lw=1.5, + color="#0072b2", + label=r"$n_2$ samples", + ) + axes_d[1].scatter( + [k2_m], + [n2_m_amse[int(len(x2_m_arr) * k2_m) - 1]], + color="#0072b2", + marker="o", + edgecolor="black", + alpha=0.5, + label=r"Min for $n_2$ sample", + ) + axes_d[1].axvline( + max_m_index1 / float(len(x1_m_arr)), + color="#d55e00", + ls="--", + alpha=0.5, + label=r"Minimization boundary for $n_1$ sample", + ) + axes_d[1].axvline( + max_m_index2 / float(len(x2_m_arr)), + color="#0072b2", + ls="--", + alpha=0.5, + label=r"Minimization boundary for $n_2$ sample", + ) + axes_d[1].legend(loc="best") if savedata and output_file_path is not None: - with open(os.path.join(output_dir+"/"+output_name+"_mom_diagn1.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_mom_diagn1.dat"), "w" + ) as f: for i in range(len(x1_m_arr[min_k1:max_k1])): - f.write(str(x1_m_arr[min_k1:max_k1][i]) + " " + str(n1_m_amse[min_k1:max_k1][i]) + "\n") - with open(os.path.join(output_dir+"/"+output_name+"_mom_diagn2.dat"), "w") as f: + f.write( + str(x1_m_arr[min_k1:max_k1][i]) + + " " + + str(n1_m_amse[min_k1:max_k1][i]) + + "\n" + ) + with open( + os.path.join(output_dir + "/" + output_name + "_mom_diagn2.dat"), "w" + ) as f: for i in range(len(x2_m_arr[min_k2:max_k2])): - f.write(str(x2_m_arr[min_k2:max_k2][i]) + " " + str(n2_m_amse[min_k2:max_k2][i]) + "\n") - with open(os.path.join(output_dir+"/"+output_name+"_mom_diagn_points.dat"), "w") as f: - f.write("Min for n1 sample: "+str(k1_m)+" "+str(n1_m_amse[int(len(x1_m_arr)*k1_m)-1])+"\n") - f.write("Min for n2 sample: "+str(k2_m)+" "+str(n2_m_amse[int(len(x2_m_arr)*k2_m)-1])+"\n") - f.write("Minimization boundary for n1 sample: "+str(max_m_index1/float(len(x1_m_arr)))+"\n") - f.write("Minimization boundary for n2 sample: "+str(max_m_index2/float(len(x2_m_arr)))+"\n") - + f.write( + str(x2_m_arr[min_k2:max_k2][i]) + + " " + + str(n2_m_amse[min_k2:max_k2][i]) + + "\n" + ) + with open( + os.path.join(output_dir + "/" + output_name + "_mom_diagn_points.dat"), "w" + ) as f: + f.write( + "Min for n1 sample: " + + str(k1_m) + + " " + + str(n1_m_amse[int(len(x1_m_arr) * k1_m) - 1]) + + "\n" + ) + f.write( + "Min for n2 sample: " + + str(k2_m) + + " " + + str(n2_m_amse[int(len(x2_m_arr) * k2_m) - 1]) + + "\n" + ) + f.write( + "Minimization boundary for n1 sample: " + + str(max_m_index1 / float(len(x1_m_arr))) + + "\n" + ) + f.write( + "Minimization boundary for n2 sample: " + + str(max_m_index2 / float(len(x2_m_arr))) + + "\n" + ) min_k1 = 2 max_k1 = len(x1_k_arr) min_k2 = 2 max_k2 = len(x2_k_arr) - n1_k_amse[np.where((n1_k_amse == np.inf) |\ - (n1_k_amse == -np.inf))] = np.nan + n1_k_amse[np.where((n1_k_amse == np.inf) | (n1_k_amse == -np.inf))] = np.nan axes_d[2].set_yscale("log") - axes_d[2].set_ylim((0.1*np.nanmin(n1_k_amse[min_k1:max_k1]), 1.0)) - axes_d[2].set_xlabel("Fraction of Bootstrap Order Statistics", - fontsize = 20) - axes_d[2].set_ylabel(r"$\langle AMSE \rangle$", fontsize = 20) - axes_d[2].set_title("Kernel-type Estimator", fontsize = 20) + axes_d[2].set_ylim((0.1 * np.nanmin(n1_k_amse[min_k1:max_k1]), 1.0)) + axes_d[2].set_xlabel("Fraction of Bootstrap Order Statistics", fontsize=20) + axes_d[2].set_ylabel(r"$\langle AMSE \rangle$", fontsize=20) + axes_d[2].set_title("Kernel-type Estimator", fontsize=20) # plot AMSE and corresponding minimum - axes_d[2].plot(x1_k_arr[min_k1:max_k1], n1_k_amse[min_k1:max_k1], - alpha = 0.5, lw = 1.5, - color = "#d55e00", label = r"$n_1$ samples") - axes_d[2].scatter([h1], [n1_k_amse[np.where(x1_k_arr == h1)][0]], color = "#d55e00", - marker = 'o', edgecolor = "black", alpha = 0.5, - label = r"Min for $n_1$ sample") + axes_d[2].plot( + x1_k_arr[min_k1:max_k1], + n1_k_amse[min_k1:max_k1], + alpha=0.5, + lw=1.5, + color="#d55e00", + label=r"$n_1$ samples", + ) + axes_d[2].scatter( + [h1], + [n1_k_amse[np.where(x1_k_arr == h1)][0]], + color="#d55e00", + marker="o", + edgecolor="black", + alpha=0.5, + label=r"Min for $n_1$ sample", + ) # plot boundary of minimization - axes_d[2].axvline(max_k_index1, color = "#d55e00", - ls = '--', alpha = 0.5, - label = r"Minimization boundary for $n_2$ sample") - axes_d[2].plot(x2_k_arr[min_k2:max_k2], n2_k_amse[min_k2:max_k2], - alpha = 0.5, lw = 1.5, - color = "#0072b2", label = r"$n_2$ samples") - axes_d[2].scatter([h2], [n2_k_amse[np.where(x2_k_arr == h2)][0]], color = "#0072b2", - marker = 'o', edgecolor = "black", alpha = 0.5, - label = r"Min for $n_2$ sample") - axes_d[2].axvline(max_k_index2, color = "#0072b2", - ls = '--', alpha = 0.5, - label = r"Minimization boundary for $n_2$ sample") - axes_d[2].legend(loc = "best") + axes_d[2].axvline( + max_k_index1, + color="#d55e00", + ls="--", + alpha=0.5, + label=r"Minimization boundary for $n_2$ sample", + ) + axes_d[2].plot( + x2_k_arr[min_k2:max_k2], + n2_k_amse[min_k2:max_k2], + alpha=0.5, + lw=1.5, + color="#0072b2", + label=r"$n_2$ samples", + ) + axes_d[2].scatter( + [h2], + [n2_k_amse[np.where(x2_k_arr == h2)][0]], + color="#0072b2", + marker="o", + edgecolor="black", + alpha=0.5, + label=r"Min for $n_2$ sample", + ) + axes_d[2].axvline( + max_k_index2, + color="#0072b2", + ls="--", + alpha=0.5, + label=r"Minimization boundary for $n_2$ sample", + ) + axes_d[2].legend(loc="best") if savedata and output_file_path is not None: - with open(os.path.join(output_dir+"/"+output_name+"_kern_diagn1.dat"), "w") as f: + with open( + os.path.join(output_dir + "/" + output_name + "_kern_diagn1.dat"), "w" + ) as f: for i in range(len(x1_k_arr[min_k1:max_k1])): - f.write(str(x1_k_arr[min_k1:max_k1][i]) + " " + str(n1_k_amse[min_k1:max_k1][i]) + "\n") - with open(os.path.join(output_dir+"/"+output_name+"_kern_diagn2.dat"), "w") as f: + f.write( + str(x1_k_arr[min_k1:max_k1][i]) + + " " + + str(n1_k_amse[min_k1:max_k1][i]) + + "\n" + ) + with open( + os.path.join(output_dir + "/" + output_name + "_kern_diagn2.dat"), "w" + ) as f: for i in range(len(x2_m_arr[min_k2:max_k2])): - f.write(str(x2_k_arr[min_k2:max_k2][i]) + " " + str(n2_k_amse[min_k2:max_k2][i]) + "\n") - with open(os.path.join(output_dir+"/"+output_name+"_kern_diagn_points.dat"), "w") as f: - f.write("Min for n1 sample: "+str(h1)+" "+str(n1_k_amse[np.where(x1_k_arr == h1)][0])+"\n") - f.write("Min for n2 sample: "+str(h2)+" "+str(n2_k_amse[np.where(x2_k_arr == h2)][0])+"\n") - f.write("Minimization boundary for n1 sample: "+str(n1_k_amse[int(max_k_index1*hsteps)-1])+"\n") - f.write("Minimization boundary for n2 sample: "+str(n2_k_amse[int(max_k_index2*hsteps)-1])+"\n") + f.write( + str(x2_k_arr[min_k2:max_k2][i]) + + " " + + str(n2_k_amse[min_k2:max_k2][i]) + + "\n" + ) + with open( + os.path.join(output_dir + "/" + output_name + "_kern_diagn_points.dat"), "w" + ) as f: + f.write( + "Min for n1 sample: " + + str(h1) + + " " + + str(n1_k_amse[np.where(x1_k_arr == h1)][0]) + + "\n" + ) + f.write( + "Min for n2 sample: " + + str(h2) + + " " + + str(n2_k_amse[np.where(x2_k_arr == h2)][0]) + + "\n" + ) + f.write( + "Minimization boundary for n1 sample: " + + str(n1_k_amse[int(max_k_index1 * hsteps) - 1]) + + "\n" + ) + f.write( + "Minimization boundary for n2 sample: " + + str(n2_k_amse[int(max_k_index2 * hsteps) - 1]) + + "\n" + ) fig_d.tight_layout() if output_file_path is not None: - diag_plots_path = os.path.dirname(os.path.realpath(output_file_path))+"/"+os.path.splitext(os.path.basename(output_file_path))[0]+"_diag.pdf" + diag_plots_path = ( + os.path.dirname(os.path.realpath(output_file_path)) + + "/" + + os.path.splitext(os.path.basename(output_file_path))[0] + + "_diag.pdf" + ) fig_d.savefig(diag_plots_path) - - return fig_d, axes_d \ No newline at end of file + + return fig_d, axes_d diff --git a/src/tailestim/estimators/result.py b/src/tailestim/estimators/result.py index ad5027b..c034428 100644 --- a/src/tailestim/estimators/result.py +++ b/src/tailestim/estimators/result.py @@ -1,4 +1,5 @@ -from typing import Dict, Any +from typing import Any, Dict + import numpy as np @@ -34,21 +35,21 @@ class TailEstimatorResult: # Mapping of keys to human-readable labels _key_labels = { - 'estimator': 'Estimator', - 'xi_star_': 'Tail index (ξ)', - 'gamma_': 'Power law exponent (γ)', - 'k_arr_': 'Order statistics', - 'xi_arr_': 'Tail index estimates', - 'k_star_': 'Optimal order statistic (k*)', - 'bootstrap_results_': 'Bootstrap Results', - 'first_bootstrap_': 'First Bootstrap', - 'second_bootstrap_': 'Second Bootstrap', - 'k_min_': 'Minimum AMSE fraction', - 'amse_': 'AMSE values', - 'max_index_': 'Maximum index', - 'x_arr_': 'Fraction of order statistics' + "estimator": "Estimator", + "xi_star_": "Tail index (ξ)", + "gamma_": "Power law exponent (γ)", + "k_arr_": "Order statistics", + "xi_arr_": "Tail index estimates", + "k_star_": "Optimal order statistic (k*)", + "bootstrap_results_": "Bootstrap Results", + "first_bootstrap_": "First Bootstrap", + "second_bootstrap_": "Second Bootstrap", + "k_min_": "Minimum AMSE fraction", + "amse_": "AMSE values", + "max_index_": "Maximum index", + "x_arr_": "Fraction of order statistics", } - + def __init__(self, initial_data: Dict[str, Any] = None, **kwargs): data = initial_data or {} data.update(kwargs) @@ -57,11 +58,11 @@ def __init__(self, initial_data: Dict[str, Any] = None, **kwargs): if isinstance(value, dict): data[key] = TailEstimatorResult(value) self.__dict__.update(data) - + def __str__(self, include_header=True) -> str: """ Return a human-readable string representation of the TailEstimatorResult. - + Parameters ---------- include_header : bool, default=True @@ -69,11 +70,11 @@ def __str__(self, include_header=True) -> str: """ output = [] from .base import BaseTailEstimator - + # Add a header for the result section, but only for the top-level object if include_header: - output.extend(["-"* 50, "Result", "-"* 50]) - + output.extend(["-" * 50, "Result", "-" * 50]) + # Order the attributes based on the keys in _key_labels ordered_keys = [k for k in self._key_labels if k in self.__dict__] # Add any keys that aren't in _key_labels at the end @@ -81,7 +82,7 @@ def __str__(self, include_header=True) -> str: for key in ordered_keys: value = self.__dict__[key] - label = self._key_labels.get(key, key.replace('_', ' ').title()) + label = self._key_labels.get(key, key.replace("_", " ").title()) # Format the value based on its type if isinstance(value, BaseTailEstimator): @@ -99,17 +100,17 @@ def __str__(self, include_header=True) -> str: elif isinstance(value, TailEstimatorResult): # For nested TailEstimatorResult objects, indent their string representation # but don't include the header - nested_str = value.__str__(include_header=False).replace('\n', '\n ') + nested_str = value.__str__(include_header=False).replace("\n", "\n ") value_str = f"\n {nested_str}" else: value_str = str(value) - + output.append(f"{label}: {value_str}") - - return '\n'.join(output) + + return "\n".join(output) def __repr__(self) -> str: """ Return a string representation of the TailEstimatorResult. """ - return f"TailEstimatorResult(estimator={self.estimator.__class__.__name__})" \ No newline at end of file + return f"TailEstimatorResult(estimator={self.estimator.__class__.__name__})" diff --git a/src/tailestim/estimators/smooth_hill.py b/src/tailestim/estimators/smooth_hill.py index 880fafa..2204c84 100644 --- a/src/tailestim/estimators/smooth_hill.py +++ b/src/tailestim/estimators/smooth_hill.py @@ -1,16 +1,20 @@ """Smooth Hill estimator implementation for tail index estimation.""" + +from typing import Any, Dict, Tuple + import numpy as np -from typing import Dict, Any, Tuple -from .result import TailEstimatorResult + from .base import BaseTailEstimator +from .result import TailEstimatorResult from .tail_methods import smooth_hill_estimator as smooth_hill_estimate + class SmoothHillEstimator(BaseTailEstimator): """Smooth Hill estimator for tail index estimation. - + This class implements the Smooth Hill estimator, which applies smoothing to the classical Hill estimator. It does not use bootstrap procedures. - + Parameters ---------- r_smooth : int, default=2 @@ -20,7 +24,7 @@ class SmoothHillEstimator(BaseTailEstimator): Additional parameters (not used by this estimator). """ - + def __init__(self, r_smooth: int = 2, **kwargs): # Smooth Hill estimator doesn't use bootstrap super().__init__(bootstrap=False, **kwargs) @@ -28,12 +32,12 @@ def __init__(self, r_smooth: int = 2, **kwargs): def _estimate(self, ordered_data: np.ndarray) -> Tuple: """Estimate the tail index using the Smooth Hill method. - + Parameters ---------- ordered_data : np.ndarray Data array in decreasing order. - + Returns ------- Tuple @@ -43,20 +47,17 @@ def _estimate(self, ordered_data: np.ndarray) -> Tuple: def get_params(self) -> Dict[str, Any]: """Get the parameters of the estimator. - + Returns ------- dict Dictionary containing the parameters of the estimator. """ - return { - "r_smooth": self.r_smooth, - **self.kwargs - } + return {"r_smooth": self.r_smooth, **self.kwargs} def get_result(self) -> TailEstimatorResult: """Get the estimated parameters. - + Attributes ---------- estimator : BaseTailEstimator @@ -72,13 +73,9 @@ def get_result(self) -> TailEstimatorResult: """ if self.results is None: raise ValueError("Model not fitted yet. Call fit() first.") - + k_arr, xi_arr = self.results - - res = { - 'estimator': self, - 'k_arr_': k_arr, - 'xi_arr_': xi_arr - } - - return TailEstimatorResult(res) \ No newline at end of file + + res = {"estimator": self, "k_arr_": k_arr, "xi_arr_": xi_arr} + + return TailEstimatorResult(res) diff --git a/src/tailestim/estimators/tail_methods.py b/src/tailestim/estimators/tail_methods.py index c797735..65df9b2 100644 --- a/src/tailestim/estimators/tail_methods.py +++ b/src/tailestim/estimators/tail_methods.py @@ -1,10 +1,12 @@ -import numpy as np -import sys import logging +import sys + +import numpy as np + logging.basicConfig(level=logging.WARNING) -def add_uniform_noise(data_sequence, p = 1, base_seed = None): +def add_uniform_noise(data_sequence, p=1, base_seed=None): """ Function to add uniform random noise to a given dataset. Uniform noise in range [-5*10^(-p), 5*10^(-p)] is added to each @@ -27,14 +29,16 @@ def add_uniform_noise(data_sequence, p = 1, base_seed = None): cur_seed = base_rng.integers(0, 1_000_000) cur_rng = np.random.default_rng(cur_seed) - noise = cur_rng.uniform(-5.*10**(-p), 5*10**(-p), size = len(data_sequence)) + noise = cur_rng.uniform(-5.0 * 10 ** (-p), 5 * 10 ** (-p), size=len(data_sequence)) randomized_data_sequence = data_sequence + noise # ensure there are no negative entries after noise is added - randomized_data_sequence = randomized_data_sequence[np.where(randomized_data_sequence > 0)] + randomized_data_sequence = randomized_data_sequence[ + np.where(randomized_data_sequence > 0) + ] return randomized_data_sequence -def get_distribution(data_sequence, number_of_bins = 30): +def get_distribution(data_sequence, number_of_bins=30): """ Function to get a log-binned distribution of a given dataset. @@ -56,21 +60,22 @@ def get_distribution(data_sequence, number_of_bins = 30): lower_bound = log(lower_bound) if lower_bound > 0 else -1 upper_bound = log(upper_bound) bins = np.logspace(lower_bound, upper_bound, number_of_bins) - + # compute the histogram using numpy - y, __ = np.histogram(data_sequence, bins = bins, density = True) + y, __ = np.histogram(data_sequence, bins=bins, density=True) # for each bin, compute its midpoint x = bins[1:] - np.diff(bins) / 2.0 # if bin is empty, drop it from the resulting list - drop_indices = [i for i,k in enumerate(y) if k == 0.0] - x = [k for i,k in enumerate(x) if i not in drop_indices] - y = [k for i,k in enumerate(y) if i not in drop_indices] + drop_indices = [i for i, k in enumerate(y) if k == 0.0] + x = [k for i, k in enumerate(x) if i not in drop_indices] + y = [k for i, k in enumerate(y) if i not in drop_indices] return x, y + def get_ccdf(degree_sequence): """ Function to get CCDF of the list of degrees. - + Args: degree_sequence: numpy array of nodes' degrees. @@ -81,8 +86,9 @@ def get_ccdf(degree_sequence): """ uniques, counts = np.unique(degree_sequence, return_counts=True) cumprob = np.cumsum(counts).astype(np.double) / (degree_sequence.size) - return uniques[::-1], (1. - cumprob)[::-1] - + return uniques[::-1], (1.0 - cumprob)[::-1] + + # ================================================ # ========== Hill Tail Index Estimation ========== # ================================================ @@ -105,142 +111,167 @@ def get_moments_estimates_1(ordered_data): logs_1 = np.log(ordered_data) logs_1_cumsum = np.cumsum(logs_1[:-1]) k_vector = np.arange(1, len(ordered_data)) - M1 = (1./k_vector)*logs_1_cumsum - logs_1[1:] + M1 = (1.0 / k_vector) * logs_1_cumsum - logs_1[1:] return M1 + def get_moments_estimates_2(ordered_data): """ Function to calculate first and second moments arrays - given an ordered data sequence. + given an ordered data sequence. Decreasing ordering is required. Args: ordered_data: numpy array of ordered data for which - the 1st (Hill estimator) and 2nd moments + the 1st (Hill estimator) and 2nd moments are calculated. Returns: M1: numpy array of 1st moments (Hill estimator) corresponding to all possible order statistics of the dataset. - M2: numpy array of 2nd moments corresponding to all + M2: numpy array of 2nd moments corresponding to all possible order statistics of the dataset. """ logs_1 = np.log(ordered_data) - logs_2 = (np.log(ordered_data))**2 + logs_2 = (np.log(ordered_data)) ** 2 logs_1_cumsum = np.cumsum(logs_1[:-1]) logs_2_cumsum = np.cumsum(logs_2[:-1]) k_vector = np.arange(1, len(ordered_data)) - M1 = (1./k_vector)*logs_1_cumsum - logs_1[1:] - M2 = (1./k_vector)*logs_2_cumsum - (2.*logs_1[1:]/k_vector)*logs_1_cumsum\ - + logs_2[1:] + M1 = (1.0 / k_vector) * logs_1_cumsum - logs_1[1:] + M2 = ( + (1.0 / k_vector) * logs_2_cumsum + - (2.0 * logs_1[1:] / k_vector) * logs_1_cumsum + + logs_2[1:] + ) return M1, M2 + def get_moments_estimates_3(ordered_data): """ - Function to calculate first, second and third moments - arrays given an ordered data sequence. + Function to calculate first, second and third moments + arrays given an ordered data sequence. Decreasing ordering is required. Args: ordered_data: numpy array of ordered data for which - the 1st (Hill estimator), 2nd and 3rd moments + the 1st (Hill estimator), 2nd and 3rd moments are calculated. Returns: M1: numpy array of 1st moments (Hill estimator) corresponding to all possible order statistics of the dataset. - M2: numpy array of 2nd moments corresponding to all + M2: numpy array of 2nd moments corresponding to all possible order statistics of the dataset. - M3: numpy array of 3rd moments corresponding to all + M3: numpy array of 3rd moments corresponding to all possible order statistics of the dataset. """ logs_1 = np.log(ordered_data) - logs_2 = (np.log(ordered_data))**2 - logs_3 = (np.log(ordered_data))**3 + logs_2 = (np.log(ordered_data)) ** 2 + logs_3 = (np.log(ordered_data)) ** 3 logs_1_cumsum = np.cumsum(logs_1[:-1]) logs_2_cumsum = np.cumsum(logs_2[:-1]) logs_3_cumsum = np.cumsum(logs_3[:-1]) k_vector = np.arange(1, len(ordered_data)) - M1 = (1./k_vector)*logs_1_cumsum - logs_1[1:] - M2 = (1./k_vector)*logs_2_cumsum - (2.*logs_1[1:]/k_vector)*logs_1_cumsum\ - + logs_2[1:] - M3 = (1./k_vector)*logs_3_cumsum - (3.*logs_1[1:]/k_vector)*logs_2_cumsum\ - + (3.*logs_2[1:]/k_vector)*logs_1_cumsum - logs_3[1:] + M1 = (1.0 / k_vector) * logs_1_cumsum - logs_1[1:] + M2 = ( + (1.0 / k_vector) * logs_2_cumsum + - (2.0 * logs_1[1:] / k_vector) * logs_1_cumsum + + logs_2[1:] + ) + M3 = ( + (1.0 / k_vector) * logs_3_cumsum + - (3.0 * logs_1[1:] / k_vector) * logs_2_cumsum + + (3.0 * logs_2[1:] / k_vector) * logs_1_cumsum + - logs_3[1:] + ) # cleaning exceptional cases - clean_indices = np.where((M2 <= 0) | (M3 == 0) | (np.abs(1.-(M1**2)/M2) < 1e-10)\ - |(np.abs(1.-(M1*M2)/M3) < 1e-10)) + clean_indices = np.where( + (M2 <= 0) + | (M3 == 0) + | (np.abs(1.0 - (M1**2) / M2) < 1e-10) + | (np.abs(1.0 - (M1 * M2) / M3) < 1e-10) + ) M1[clean_indices] = np.nan M2[clean_indices] = np.nan M3[clean_indices] = np.nan return M1, M2, M3 -def hill_dbs(ordered_data, t_bootstrap = 0.5, - r_bootstrap = 500, eps_stop = 1.0, - verbose = False, diagn_plots = False, base_seed=None, - max_resample = 50): + +def hill_dbs( + ordered_data, + t_bootstrap=0.5, + r_bootstrap=500, + eps_stop=1.0, + verbose=False, + diagn_plots=False, + base_seed=None, + max_resample=50, +): """ - Function to perform double-bootstrap procedure for - Hill estimator. - - Args: - ordered_data: numpy array for which double-bootstrap - is performed. Decreasing ordering is required. - t_bootstrap: parameter controlling the size of the 2nd - bootstrap. Defined from n2 = n*(t_bootstrap). - r_bootstrap: number of bootstrap resamplings for the 1st and 2nd - bootstraps. - eps_stop: parameter controlling range of AMSE minimization. - Defined as the fraction of order statistics to consider - during the AMSE minimization step. - verbose: flag controlling bootstrap verbosity. - diagn_plots: flag to switch on/off generation of AMSE diagnostic - plots. - base_seed: base random seed for reproducibility of bootstrap (default is None). - max_resample: maximum number of resampling attempts when AMSE - false minimum is detected (k2 > k1). Raises - ConvergenceError if exceeded. Default is 50. - - Returns: - k_star: number of order statistics optimal for estimation - according to the double-bootstrap procedure. - x1_arr: array of fractions of order statistics used for the - 1st bootstrap sample. - n1_amse: array of AMSE values produced by the 1st bootstrap - sample. - k1_min: value of fraction of order statistics corresponding - to the minimum of AMSE for the 1st bootstrap sample. - max_index1: index of the 1st bootstrap sample's order statistics - array corresponding to the minimization boundary set - by eps_stop parameter. - x2_arr: array of fractions of order statistics used for the - 2nd bootstrap sample. - n2_amse: array of AMSE values produced by the 2nd bootstrap - sample. - k2_min: value of fraction of order statistics corresponding - to the minimum of AMSE for the 2nd bootstrap sample. - max_index2: index of the 2nd bootstrap sample's order statistics - array corresponding to the minimization boundary set - by eps_stop parameter. + Function to perform double-bootstrap procedure for + Hill estimator. + + Args: + ordered_data: numpy array for which double-bootstrap + is performed. Decreasing ordering is required. + t_bootstrap: parameter controlling the size of the 2nd + bootstrap. Defined from n2 = n*(t_bootstrap). + r_bootstrap: number of bootstrap resamplings for the 1st and 2nd + bootstraps. + eps_stop: parameter controlling range of AMSE minimization. + Defined as the fraction of order statistics to consider + during the AMSE minimization step. + verbose: flag controlling bootstrap verbosity. + diagn_plots: flag to switch on/off generation of AMSE diagnostic + plots. + base_seed: base random seed for reproducibility of bootstrap (default is None). + max_resample: maximum number of resampling attempts when AMSE + false minimum is detected (k2 > k1). Raises + ConvergenceError if exceeded. Default is 50. + + Returns: + k_star: number of order statistics optimal for estimation + according to the double-bootstrap procedure. + x1_arr: array of fractions of order statistics used for the + 1st bootstrap sample. + n1_amse: array of AMSE values produced by the 1st bootstrap + sample. + k1_min: value of fraction of order statistics corresponding + to the minimum of AMSE for the 1st bootstrap sample. + max_index1: index of the 1st bootstrap sample's order statistics + array corresponding to the minimization boundary set + by eps_stop parameter. + x2_arr: array of fractions of order statistics used for the + 2nd bootstrap sample. + n2_amse: array of AMSE values produced by the 2nd bootstrap + sample. + k2_min: value of fraction of order statistics corresponding + to the minimum of AMSE for the 2nd bootstrap sample. + max_index2: index of the 2nd bootstrap sample's order statistics + array corresponding to the minimization boundary set + by eps_stop parameter. """ if verbose: logging.debug("Performing Hill double-bootstrap...") n = len(ordered_data) - eps_bootstrap = 0.5*(1+np.log(int(t_bootstrap*n))/np.log(n)) + eps_bootstrap = 0.5 * (1 + np.log(int(t_bootstrap * n)) / np.log(n)) n1 = int(n**eps_bootstrap) - samples_n1 = np.zeros(n1-1) - good_counts1 = np.zeros(n1-1) + samples_n1 = np.zeros(n1 - 1) + good_counts1 = np.zeros(n1 - 1) k1 = None k2 = None min_index1 = 1 min_index2 = 1 - base_rng = np.random.default_rng(seed=base_seed) # Accept random seed for reproducibility. Default seed is None. + base_rng = np.random.default_rng( + seed=base_seed + ) # Accept random seed for reproducibility. Default seed is None. resample_count = 0 - while k2 == None: + while k2 is None: if resample_count >= max_resample: raise RuntimeError( f"Hill double-bootstrap failed to converge after " @@ -248,73 +279,84 @@ def hill_dbs(ordered_data, t_bootstrap = 0.5, f"Consider increasing max_resample or adjusting bootstrap parameters." ) # first bootstrap with n1 sample size - for i in range(r_bootstrap): + for _i in range(r_bootstrap): bs1_seed = base_rng.integers(0, 1_000_000) bs1_rng = np.random.default_rng(bs1_seed) - sample = bs1_rng.choice(ordered_data, n1, replace = True) + sample = bs1_rng.choice(ordered_data, n1, replace=True) sample[::-1].sort() M1, M2 = get_moments_estimates_2(sample) - current_amse1 = (M2 - 2.*(M1)**2)**2 + current_amse1 = (M2 - 2.0 * (M1) ** 2) ** 2 samples_n1 += current_amse1 good_counts1[np.where(current_amse1 != np.nan)] += 1 averaged_delta = samples_n1 / good_counts1 - - max_index1 = (np.abs(np.linspace(1./n1, 1.0, n1) - eps_stop)).argmin() - k1 = np.nanargmin(averaged_delta[min_index1:max_index1]) + 1 + min_index1 #take care of indexing + + max_index1 = (np.abs(np.linspace(1.0 / n1, 1.0, n1) - eps_stop)).argmin() + k1 = ( + np.nanargmin(averaged_delta[min_index1:max_index1]) + 1 + min_index1 + ) # take care of indexing if diagn_plots: n1_amse = averaged_delta - x1_arr = np.linspace(1./n1, 1.0, n1) - + x1_arr = np.linspace(1.0 / n1, 1.0, n1) + # second bootstrap with n2 sample size - n2 = int(n1*n1/float(n)) - samples_n2 = np.zeros(n2-1) - good_counts2 = np.zeros(n2-1) - - for i in range(r_bootstrap): + n2 = int(n1 * n1 / float(n)) + samples_n2 = np.zeros(n2 - 1) + good_counts2 = np.zeros(n2 - 1) + + for _i in range(r_bootstrap): bs2_seed = base_rng.integers(0, 1_000_000) bs2_rng = np.random.default_rng(bs2_seed) - sample = bs2_rng.choice(ordered_data, n2, replace = True) + sample = bs2_rng.choice(ordered_data, n2, replace=True) sample[::-1].sort() M1, M2 = get_moments_estimates_2(sample) - current_amse2 = (M2 - 2.*(M1**2))**2 + current_amse2 = (M2 - 2.0 * (M1**2)) ** 2 samples_n2 += current_amse2 good_counts2[np.where(current_amse2 != np.nan)] += 1 - max_index2 = (np.abs(np.linspace(1./n2, 1.0, n2) - eps_stop)).argmin() + max_index2 = (np.abs(np.linspace(1.0 / n2, 1.0, n2) - eps_stop)).argmin() averaged_delta = samples_n2 / good_counts2 - - - k2 = np.nanargmin(averaged_delta[min_index2:max_index2]) + 1 + min_index2 #take care of indexing + + k2 = ( + np.nanargmin(averaged_delta[min_index2:max_index2]) + 1 + min_index2 + ) # take care of indexing if diagn_plots: n2_amse = averaged_delta - x2_arr = np.linspace(1./n2, 1.0, n2) + x2_arr = np.linspace(1.0 / n2, 1.0, n2) if k2 > k1: resample_count += 1 - logging.warning("Warning (Hill): k2 > k1, AMSE false minimum suspected, resampling... (attempt %d/%d)", resample_count, max_resample) + logging.warning( + "Warning (Hill): k2 > k1, AMSE false minimum suspected, resampling... (attempt %d/%d)", + resample_count, + max_resample, + ) # move left AMSE boundary to avoid numerical issues - min_index1 = min_index1 + int(0.005*n) - min_index2 = min_index2 + int(0.005*n) + min_index1 = min_index1 + int(0.005 * n) + min_index2 = min_index2 + int(0.005 * n) k2 = None - - ''' + + """ # this constant is provided in the Danielsson's paper # use instead of rho below if needed rho = (np.log(k1)/(2.*np.log(n1) - np.log(k1)))\ **(2.*(np.log(n1) - np.log(k1))/(np.log(n1))) - ''' - + """ + # this constant is provided in Qi's paper - rho = (1. - (2*(np.log(k1) - np.log(n1))/(np.log(k1))))**(np.log(k1)/np.log(n1) - 1.) - - k_star = (k1*k1/float(k2)) * rho + rho = (1.0 - (2 * (np.log(k1) - np.log(n1)) / (np.log(k1)))) ** ( + np.log(k1) / np.log(n1) - 1.0 + ) + + k_star = (k1 * k1 / float(k2)) * rho k_star = int(np.round(k_star)) - + # enforce k_star to pick 2nd value (rare cases of extreme cutoffs) if k_star == 0: k_star = 2 if int(k_star) >= len(ordered_data): - logging.warning("WARNING: estimated threshold k is larger than the size of data") - k_star = len(ordered_data)-1 + logging.warning( + "WARNING: estimated threshold k is larger than the size of data" + ) + k_star = len(ordered_data) - 1 if verbose: logging.info("--- Hill double-bootstrap information ---") logging.info("Size of the 1st bootstrap sample n1:", n1) @@ -326,13 +368,30 @@ def hill_dbs(ordered_data, t_bootstrap = 0.5, logging.info("-----------------------------------------") if not diagn_plots: x1_arr, x2_arr, n1_amse, n2_amse = None, None, None, None - return k_star, x1_arr, n1_amse, k1/float(n1), max_index1, x2_arr, n2_amse, k2/float(n2), max_index2 - -def hill_estimator(ordered_data, - bootstrap = True, t_bootstrap = 0.5, - r_bootstrap = 500, verbose = False, - diagn_plots = False, eps_stop = 0.99, base_seed = None, - max_resample = 50): + return ( + k_star, + x1_arr, + n1_amse, + k1 / float(n1), + max_index1, + x2_arr, + n2_amse, + k2 / float(n2), + max_index2, + ) + + +def hill_estimator( + ordered_data, + bootstrap=True, + t_bootstrap=0.5, + r_bootstrap=500, + verbose=False, + diagn_plots=False, + eps_stop=0.99, + base_seed=None, + max_resample=50, +): """ Function to calculate Hill estimator for a given dataset. If bootstrap flag is True, double-bootstrap procedure @@ -350,7 +409,7 @@ def hill_estimator(ordered_data, eps_stop: parameter controlling range of AMSE minimization. Defined as the fraction of order statistics to consider during the AMSE minimization step. - verbose: flag controlling bootstrap verbosity. + verbose: flag controlling bootstrap verbosity. diagn_plots: flag to switch on/off generation of AMSE diagnostic plots. base_seed: base random seed for reproducibility of bootstrap (default is None). @@ -374,37 +433,58 @@ def hill_estimator(ordered_data, k_arr = np.arange(1, len(ordered_data)) xi_arr = get_moments_estimates_1(ordered_data) if bootstrap: - results = hill_dbs(ordered_data, - t_bootstrap = t_bootstrap, - r_bootstrap = r_bootstrap, - verbose = verbose, - diagn_plots = diagn_plots, - eps_stop = eps_stop, - base_seed = base_seed, - max_resample = max_resample) - k_star, x1_arr, n1_amse, k1, max_index1, x2_arr, n2_amse, k2, max_index2 = results - while k_star == None: + results = hill_dbs( + ordered_data, + t_bootstrap=t_bootstrap, + r_bootstrap=r_bootstrap, + verbose=verbose, + diagn_plots=diagn_plots, + eps_stop=eps_stop, + base_seed=base_seed, + max_resample=max_resample, + ) + k_star, x1_arr, n1_amse, k1, max_index1, x2_arr, n2_amse, k2, max_index2 = ( + results + ) + while k_star is None: logging.debug("Resampling...") - results = hill_dbs(ordered_data, - t_bootstrap = t_bootstrap, - r_bootstrap = r_bootstrap, - verbose = verbose, - diagn_plots = diagn_plots, - eps_stop = eps_stop, - base_seed = base_seed, - max_resample = max_resample) - k_star, x1_arr, n1_amse, k1, max_index1, x2_arr, n2_amse, k2, max_index2 = results - xi_star = xi_arr[k_star-1] - logging.info("Adjusted Hill estimated gamma:", 1 + 1./xi_star) + results = hill_dbs( + ordered_data, + t_bootstrap=t_bootstrap, + r_bootstrap=r_bootstrap, + verbose=verbose, + diagn_plots=diagn_plots, + eps_stop=eps_stop, + base_seed=base_seed, + max_resample=max_resample, + ) + k_star, x1_arr, n1_amse, k1, max_index1, x2_arr, n2_amse, k2, max_index2 = ( + results + ) + xi_star = xi_arr[k_star - 1] + logging.info("Adjusted Hill estimated gamma:", 1 + 1.0 / xi_star) else: k_star, xi_star = None, None - x1_arr, n1_amse, k1, max_index1 = 4*[None] - x2_arr, n2_amse, k2, max_index2 = 4*[None] - results = [k_arr, xi_arr, k_star, xi_star, x1_arr, n1_amse, k1, max_index1,\ - x2_arr, n2_amse, k2, max_index2] + x1_arr, n1_amse, k1, max_index1 = 4 * [None] + x2_arr, n2_amse, k2, max_index2 = 4 * [None] + results = [ + k_arr, + xi_arr, + k_star, + xi_star, + x1_arr, + n1_amse, + k1, + max_index1, + x2_arr, + n2_amse, + k2, + max_index2, + ] return results -def smooth_hill_estimator(ordered_data, r_smooth = 2): + +def smooth_hill_estimator(ordered_data, r_smooth=2): """ Function to calculate smooth Hill estimator for a given ordered dataset. @@ -417,31 +497,33 @@ def smooth_hill_estimator(ordered_data, r_smooth = 2): value such as 2 or 3. Returns: k_arr: numpy array of order statistics based on the data provided. - xi_arr: numpy array of tail index estimates corresponding to + xi_arr: numpy array of tail index estimates corresponding to the order statistics array k_arr. """ n = len(ordered_data) M1 = get_moments_estimates_1(ordered_data) - xi_arr = np.zeros(int(np.floor(float(n)/r_smooth))) - k_arr = np.arange(1, int(np.floor(float(n)/r_smooth))+1) + xi_arr = np.zeros(int(np.floor(float(n) / r_smooth))) + k_arr = np.arange(1, int(np.floor(float(n) / r_smooth)) + 1) xi_arr[0] = M1[0] - bin_lengths = np.array([1.]+[float((r_smooth-1)*k) for k in k_arr[:-1]]) + bin_lengths = np.array([1.0] + [float((r_smooth - 1) * k) for k in k_arr[:-1]]) cum_sum = 0.0 - for i in range(1, r_smooth*int(np.floor(float(n)/r_smooth))-1): + for i in range(1, r_smooth * int(np.floor(float(n) / r_smooth)) - 1): k = i cum_sum += M1[k] - if (k+1) % (r_smooth) == 0: - xi_arr[int(k+1)//int(r_smooth)] = cum_sum - cum_sum -= M1[int(k+1)//int(r_smooth)] - xi_arr = xi_arr/bin_lengths + if (k + 1) % (r_smooth) == 0: + xi_arr[int(k + 1) // int(r_smooth)] = cum_sum + cum_sum -= M1[int(k + 1) // int(r_smooth)] + xi_arr = xi_arr / bin_lengths return k_arr, xi_arr + # =================================================== # ========== Moments Tail Index Estimation ========== # =================================================== + def moments_dbs_prefactor(xi_n, n1, k1): - """ + """ Function to calculate pre-factor used in moments double-bootstrap procedure. @@ -457,65 +539,74 @@ def moments_dbs_prefactor(xi_n, n1, k1): prefactor: constant used in estimation of the optimal stopping order statistic for moments estimator. """ + def V_sq(xi_n): if xi_n >= 0: - V = 1. + (xi_n)**2 + V = 1.0 + (xi_n) ** 2 return V else: - a = (1.-xi_n)**2 - b = (1-2*xi_n)*(6*((xi_n)**2)-xi_n+1) - c = (1.-3*xi_n)*(1-4*xi_n) - V = a*b/c + a = (1.0 - xi_n) ** 2 + b = (1 - 2 * xi_n) * (6 * ((xi_n) ** 2) - xi_n + 1) + c = (1.0 - 3 * xi_n) * (1 - 4 * xi_n) + V = a * b / c return V def V_bar_sq(xi_n): if xi_n >= 0: - V = 0.25*(1+(xi_n)**2) + V = 0.25 * (1 + (xi_n) ** 2) return V else: - a = 0.25*((1-xi_n)**2) - b = 1-8*xi_n+48*(xi_n**2)-154*(xi_n**3) - c = 263*(xi_n**4)-222*(xi_n**5)+72*(xi_n**6) - d = (1.-2*xi_n)*(1-3*xi_n)*(1-4*xi_n) - e = (1.-5*xi_n)*(1-6*xi_n) - V = a*(b+c)/(d*e) + a = 0.25 * ((1 - xi_n) ** 2) + b = 1 - 8 * xi_n + 48 * (xi_n**2) - 154 * (xi_n**3) + c = 263 * (xi_n**4) - 222 * (xi_n**5) + 72 * (xi_n**6) + d = (1.0 - 2 * xi_n) * (1 - 3 * xi_n) * (1 - 4 * xi_n) + e = (1.0 - 5 * xi_n) * (1 - 6 * xi_n) + V = a * (b + c) / (d * e) return V - + def b(xi_n, rho): if xi_n < rho: - a1 = (1.-xi_n)*(1-2*xi_n) - a2 = (1.-rho-xi_n)*(1.-rho-2*xi_n) - return a1/a2 + a1 = (1.0 - xi_n) * (1 - 2 * xi_n) + a2 = (1.0 - rho - xi_n) * (1.0 - rho - 2 * xi_n) + return a1 / a2 elif xi_n >= rho and xi_n < 0: - return 1./(1-xi_n) + return 1.0 / (1 - xi_n) else: - b = (xi_n/(rho*(1.-rho))) + (1./((1-rho)**2)) + b = (xi_n / (rho * (1.0 - rho))) + (1.0 / ((1 - rho) ** 2)) return b def b_bar(xi_n, rho): if xi_n < rho: - a1 = 0.5*(-rho*(1-xi_n)**2) - a2 = (1.-xi_n-rho)*(1-2*xi_n-rho)*(1-3*xi_n-rho) - return a1/a2 + a1 = 0.5 * (-rho * (1 - xi_n) ** 2) + a2 = (1.0 - xi_n - rho) * (1 - 2 * xi_n - rho) * (1 - 3 * xi_n - rho) + return a1 / a2 elif xi_n >= rho and xi_n < 0: - a1 = 1-2*xi_n-np.sqrt((1-xi_n)*(1-2.*xi_n)) - a2 = (1.-xi_n)*(1-2*xi_n) - return a1/a2 + a1 = 1 - 2 * xi_n - np.sqrt((1 - xi_n) * (1 - 2.0 * xi_n)) + a2 = (1.0 - xi_n) * (1 - 2 * xi_n) + return a1 / a2 else: - b = (-1.)*((rho + xi_n*(1-rho))/(2*(1-rho)**3)) + b = (-1.0) * ((rho + xi_n * (1 - rho)) / (2 * (1 - rho) ** 3)) return b - rho = np.log(k1)/(2*np.log(k1) - 2.*np.log(n1)) - a = (V_sq(xi_n)) * (b_bar(xi_n, rho)**2) - b = V_bar_sq(xi_n) * (b(xi_n, rho)**2) - prefactor = (a/b)**(1./(1. - 2*rho)) + rho = np.log(k1) / (2 * np.log(k1) - 2.0 * np.log(n1)) + a = (V_sq(xi_n)) * (b_bar(xi_n, rho) ** 2) + b = V_bar_sq(xi_n) * (b(xi_n, rho) ** 2) + prefactor = (a / b) ** (1.0 / (1.0 - 2 * rho)) return prefactor -def moments_dbs(ordered_data, xi_n, t_bootstrap = 0.5, - r_bootstrap = 500, eps_stop = 1.0, - verbose = False, diagn_plots = False, base_seed=None): + +def moments_dbs( + ordered_data, + xi_n, + t_bootstrap=0.5, + r_bootstrap=500, + eps_stop=1.0, + verbose=False, + diagn_plots=False, + base_seed=None, +): """ - Function to perform double-bootstrap procedure for + Function to perform double-bootstrap procedure for moments estimator. Args: @@ -530,11 +621,11 @@ def moments_dbs(ordered_data, xi_n, t_bootstrap = 0.5, eps_stop: parameter controlling range of AMSE minimization. Defined as the fraction of order statistics to consider during the AMSE minimization step. - verbose: flag controlling bootstrap verbosity. + verbose: flag controlling bootstrap verbosity. diagn_plots: flag to switch on/off generation of AMSE diagnostic plots. base_seed: base random seed for reproducibility of bootstrap (default is None). - + Returns: k_star: number of order statistics optimal for estimation @@ -561,64 +652,69 @@ def moments_dbs(ordered_data, xi_n, t_bootstrap = 0.5, if verbose: logging.debug("Performing moments double-bootstrap...") n = len(ordered_data) - eps_bootstrap = 0.5*(1+np.log(int(t_bootstrap*n))/np.log(n)) + eps_bootstrap = 0.5 * (1 + np.log(int(t_bootstrap * n)) / np.log(n)) - base_rng = np.random.default_rng(seed=base_seed) # Accept random seed for reproducibility. Default seed is None. + base_rng = np.random.default_rng( + seed=base_seed + ) # Accept random seed for reproducibility. Default seed is None. # first bootstrap with n1 sample size n1 = int(n**eps_bootstrap) - samples_n1 = np.zeros(n1-1) - good_counts1 = np.zeros(n1-1) - for i in range(r_bootstrap): + samples_n1 = np.zeros(n1 - 1) + good_counts1 = np.zeros(n1 - 1) + for _i in range(r_bootstrap): bs1_seed = base_rng.integers(0, 1_000_000) bs1_rng = np.random.default_rng(bs1_seed) - sample = bs1_rng.choice(ordered_data, n1, replace = True) + sample = bs1_rng.choice(ordered_data, n1, replace=True) sample[::-1].sort() M1, M2, M3 = get_moments_estimates_3(sample) - xi_2 = M1 + 1. - 0.5*((1. - (M1*M1)/M2))**(-1.) - xi_3 = np.sqrt(0.5*M2) + 1. - (2./3.)*(1. / (1. - M1*M2/M3)) - samples_n1 += (xi_2 - xi_3)**2 - good_counts1[np.where((xi_2 - xi_3)**2 != np.nan)] += 1 - max_index1 = (np.abs(np.linspace(1./n1, 1.0, n1) - eps_stop)).argmin() + xi_2 = M1 + 1.0 - 0.5 * (1.0 - (M1 * M1) / M2) ** (-1.0) + xi_3 = np.sqrt(0.5 * M2) + 1.0 - (2.0 / 3.0) * (1.0 / (1.0 - M1 * M2 / M3)) + samples_n1 += (xi_2 - xi_3) ** 2 + good_counts1[np.where((xi_2 - xi_3) ** 2 != np.nan)] += 1 + max_index1 = (np.abs(np.linspace(1.0 / n1, 1.0, n1) - eps_stop)).argmin() averaged_delta = samples_n1 / good_counts1 - k1 = np.nanargmin(averaged_delta[:max_index1]) + 1 #take care of indexing + k1 = np.nanargmin(averaged_delta[:max_index1]) + 1 # take care of indexing if diagn_plots: n1_amse = averaged_delta - x1_arr = np.linspace(1./n1, 1.0, n1) - - - #r second bootstrap with n2 sample size - n2 = int(n1*n1/float(n)) - samples_n2 = np.zeros(n2-1) - good_counts2 = np.zeros(n2-1) - for i in range(r_bootstrap): + x1_arr = np.linspace(1.0 / n1, 1.0, n1) + + # r second bootstrap with n2 sample size + n2 = int(n1 * n1 / float(n)) + samples_n2 = np.zeros(n2 - 1) + good_counts2 = np.zeros(n2 - 1) + for _i in range(r_bootstrap): bs2_seed = base_rng.integers(0, 1_000_000) bs2_rng = np.random.default_rng(bs2_seed) - sample = bs2_rng.choice(ordered_data, n2, replace = True) + sample = bs2_rng.choice(ordered_data, n2, replace=True) sample[::-1].sort() M1, M2, M3 = get_moments_estimates_3(sample) - xi_2 = M1 + 1. - 0.5*(1. - (M1*M1)/M2)**(-1.) - xi_3 = np.sqrt(0.5*M2) + 1. - (2./3.)*(1. / (1. - M1*M2/M3)) - samples_n2 += (xi_2 - xi_3)**2 - good_counts2[np.where((xi_2 - xi_3)**2 != np.nan)] += 1 - max_index2 = (np.abs(np.linspace(1./n2, 1.0, n2) - eps_stop)).argmin() + xi_2 = M1 + 1.0 - 0.5 * (1.0 - (M1 * M1) / M2) ** (-1.0) + xi_3 = np.sqrt(0.5 * M2) + 1.0 - (2.0 / 3.0) * (1.0 / (1.0 - M1 * M2 / M3)) + samples_n2 += (xi_2 - xi_3) ** 2 + good_counts2[np.where((xi_2 - xi_3) ** 2 != np.nan)] += 1 + max_index2 = (np.abs(np.linspace(1.0 / n2, 1.0, n2) - eps_stop)).argmin() averaged_delta = samples_n2 / good_counts2 - k2 = np.nanargmin(averaged_delta[:max_index2]) + 1 #take care of indexing + k2 = np.nanargmin(averaged_delta[:max_index2]) + 1 # take care of indexing if diagn_plots: n2_amse = averaged_delta - x2_arr = np.linspace(1./n2, 1.0, n2) - + x2_arr = np.linspace(1.0 / n2, 1.0, n2) + if k2 > k1: - logging.warning("WARNING(moments): estimated k2 is greater than k1! Re-doing bootstrap...") - return 9*[None] - - #calculate estimated optimal stopping k + logging.warning( + "WARNING(moments): estimated k2 is greater than k1! Re-doing bootstrap..." + ) + return 9 * [None] + + # calculate estimated optimal stopping k prefactor = moments_dbs_prefactor(xi_n, n1, k1) - k_star = int((k1*k1/float(k2)) * prefactor) + k_star = int((k1 * k1 / float(k2)) * prefactor) if int(k_star) >= len(ordered_data): - logging.warning("WARNING: estimated threshold k is larger than the size of data") - k_star = len(ordered_data)-1 + logging.warning( + "WARNING: estimated threshold k is larger than the size of data" + ) + k_star = len(ordered_data) - 1 if verbose: logging.info("--- Moments double-bootstrap information ---") logging.info("Size of the 1st bootstrap sample n1:", n1) @@ -630,12 +726,29 @@ def moments_dbs(ordered_data, xi_n, t_bootstrap = 0.5, logging.info("--------------------------------------------") if not diagn_plots: x1_arr, x2_arr, n1_amse, n2_amse = None, None, None, None - return k_star, x1_arr, n1_amse, k1/float(n1), max_index1, x2_arr, n2_amse, k2/float(n2), max_index2 - -def moments_estimator(ordered_data, - bootstrap = True, t_bootstrap = 0.5, - r_bootstrap = 500, verbose = False, - diagn_plots = False, eps_stop = 0.99, base_seed = None): + return ( + k_star, + x1_arr, + n1_amse, + k1 / float(n1), + max_index1, + x2_arr, + n2_amse, + k2 / float(n2), + max_index2, + ) + + +def moments_estimator( + ordered_data, + bootstrap=True, + t_bootstrap=0.5, + r_bootstrap=500, + verbose=False, + diagn_plots=False, + eps_stop=0.99, + base_seed=None, +): """ Function to calculate moments estimator for a given dataset. If bootstrap flag is True, double-bootstrap procedure @@ -653,7 +766,7 @@ def moments_estimator(ordered_data, eps_stop: parameter controlling range of AMSE minimization. Defined as the fraction of order statistics to consider during the AMSE minimization step. - verbose: flag controlling bootstrap verbosity. + verbose: flag controlling bootstrap verbosity. diagn_plots: flag to switch on/off generation of AMSE diagnostic plots. base_seed: base random seed for reproducibility of bootstrap (default is None). @@ -672,46 +785,68 @@ def moments_estimator(ordered_data, by eps_stop parameter; and the same characteristics for the 2nd bootstrap sample. """ - n = len(ordered_data) + n = len(ordered_data) M1, M2 = get_moments_estimates_2(ordered_data) - xi_arr = M1 + 1. - 0.5*(1. - (M1*M1)/M2)**(-1) + xi_arr = M1 + 1.0 - 0.5 * (1.0 - (M1 * M1) / M2) ** (-1) k_arr = np.arange(1, len(ordered_data)) if bootstrap: - xi_n = xi_arr[int(np.floor(n**0.5))-1] - results = moments_dbs(ordered_data, xi_n, - t_bootstrap = t_bootstrap, - r_bootstrap = r_bootstrap, - verbose = verbose, - diagn_plots = diagn_plots, - eps_stop = eps_stop, - base_seed = base_seed) - while results[0] == None: + xi_n = xi_arr[int(np.floor(n**0.5)) - 1] + results = moments_dbs( + ordered_data, + xi_n, + t_bootstrap=t_bootstrap, + r_bootstrap=r_bootstrap, + verbose=verbose, + diagn_plots=diagn_plots, + eps_stop=eps_stop, + base_seed=base_seed, + ) + while results[0] is None: logging.debug("Resampling...") - results = moments_dbs(ordered_data, xi_n, - t_bootstrap = t_bootstrap, - r_bootstrap = r_bootstrap, - verbose = verbose, - diagn_plots = diagn_plots, - eps_stop = eps_stop, - base_seed = base_seed) - k_star, x1_arr, n1_amse, k1, max_index1, x2_arr, n2_amse, k2, max_index2 = results - xi_star = xi_arr[k_star-1] + results = moments_dbs( + ordered_data, + xi_n, + t_bootstrap=t_bootstrap, + r_bootstrap=r_bootstrap, + verbose=verbose, + diagn_plots=diagn_plots, + eps_stop=eps_stop, + base_seed=base_seed, + ) + k_star, x1_arr, n1_amse, k1, max_index1, x2_arr, n2_amse, k2, max_index2 = ( + results + ) + xi_star = xi_arr[k_star - 1] if xi_star <= 0: logging.info("Moments estimated gamma: infinity (xi <= 0).") else: - logging.info("Moments estimated gamma:", 1 + 1./xi_star) + logging.info("Moments estimated gamma:", 1 + 1.0 / xi_star) else: k_star, xi_star = None, None - x1_arr, n1_amse, k1, max_index1 = 4*[None] - x2_arr, n2_amse, k2, max_index2 = 4*[None] - results = [k_arr, xi_arr, k_star, xi_star, x1_arr, n1_amse, k1, max_index1,\ - x2_arr, n2_amse, k2, max_index2] + x1_arr, n1_amse, k1, max_index1 = 4 * [None] + x2_arr, n2_amse, k2, max_index2 = 4 * [None] + results = [ + k_arr, + xi_arr, + k_star, + xi_star, + x1_arr, + n1_amse, + k1, + max_index1, + x2_arr, + n2_amse, + k2, + max_index2, + ] return results + # ======================================================= # ========== Kernel-type Tail Index Estimation ========== # ======================================================= + def get_biweight_kernel_estimates(ordered_data, hsteps, alpha): """ Function to calculate biweight kernel-type estimates for tail index. @@ -737,33 +872,39 @@ def get_biweight_kernel_estimates(ordered_data, hsteps, alpha): n = len(ordered_data) logs = np.log(ordered_data) differences = logs[:-1] - logs[1:] - i_arr = np.arange(1, n)/float(n) + i_arr = np.arange(1, n) / float(n) i3_arr = i_arr**3 i5_arr = i_arr**5 i_alpha_arr = i_arr**alpha - i_alpha2_arr = i_arr**(2.+alpha) - i_alpha4_arr = i_arr**(4.+alpha) - t1 = np.cumsum(i_arr*differences) - t2 = np.cumsum(i3_arr*differences) - t3 = np.cumsum(i5_arr*differences) - t4 = np.cumsum(i_alpha_arr*differences) - t5 = np.cumsum(i_alpha2_arr*differences) - t6 = np.cumsum(i_alpha4_arr*differences) - h_arr = np.logspace(np.log10(1./n), np.log10(1.0), hsteps) - max_i_vector = (np.floor((n*h_arr))-2.).astype(int) - gamma_pos = (15./(8*h_arr))*t1[max_i_vector]\ - - (15./(4*(h_arr**3)))*t2[max_i_vector]\ - + (15./(8*(h_arr**5)))*t3[max_i_vector] - - q1 = (15./(8*h_arr))*t4[max_i_vector]\ - + (15./(8*(h_arr**5)))*t6[max_i_vector]\ - - (15./(4*(h_arr**3)))*t5[max_i_vector] - - q2 = (15.*(1+alpha)/(8*h_arr))*t4[max_i_vector]\ - + (15.*(5+alpha)/(8*(h_arr**5)))*t6[max_i_vector]\ - - (15.*(3+alpha)/(4*(h_arr**3)))*t5[max_i_vector] - - xi_arr = gamma_pos -1. + q2/q1 + i_alpha2_arr = i_arr ** (2.0 + alpha) + i_alpha4_arr = i_arr ** (4.0 + alpha) + t1 = np.cumsum(i_arr * differences) + t2 = np.cumsum(i3_arr * differences) + t3 = np.cumsum(i5_arr * differences) + t4 = np.cumsum(i_alpha_arr * differences) + t5 = np.cumsum(i_alpha2_arr * differences) + t6 = np.cumsum(i_alpha4_arr * differences) + h_arr = np.logspace(np.log10(1.0 / n), np.log10(1.0), hsteps) + max_i_vector = (np.floor(n * h_arr) - 2.0).astype(int) + gamma_pos = ( + (15.0 / (8 * h_arr)) * t1[max_i_vector] + - (15.0 / (4 * (h_arr**3))) * t2[max_i_vector] + + (15.0 / (8 * (h_arr**5))) * t3[max_i_vector] + ) + + q1 = ( + (15.0 / (8 * h_arr)) * t4[max_i_vector] + + (15.0 / (8 * (h_arr**5))) * t6[max_i_vector] + - (15.0 / (4 * (h_arr**3))) * t5[max_i_vector] + ) + + q2 = ( + (15.0 * (1 + alpha) / (8 * h_arr)) * t4[max_i_vector] + + (15.0 * (5 + alpha) / (8 * (h_arr**5))) * t6[max_i_vector] + - (15.0 * (3 + alpha) / (4 * (h_arr**3))) * t5[max_i_vector] + ) + + xi_arr = gamma_pos - 1.0 + q2 / q1 return h_arr, xi_arr @@ -792,48 +933,63 @@ def get_triweight_kernel_estimates(ordered_data, hsteps, alpha): n = len(ordered_data) logs = np.log(ordered_data) differences = logs[:-1] - logs[1:] - i_arr = np.arange(1, n)/float(n) + i_arr = np.arange(1, n) / float(n) i3_arr = i_arr**3 i5_arr = i_arr**5 i7_arr = i_arr**7 i_alpha_arr = i_arr**alpha - i_alpha2_arr = i_arr**(2.+alpha) - i_alpha4_arr = i_arr**(4.+alpha) - i_alpha6_arr = i_arr**(6.+alpha) - t1 = np.cumsum(i_arr*differences) - t2 = np.cumsum(i3_arr*differences) - t3 = np.cumsum(i5_arr*differences) - t4 = np.cumsum(i7_arr*differences) - t5 = np.cumsum(i_alpha_arr*differences) - t6 = np.cumsum(i_alpha2_arr*differences) - t7 = np.cumsum(i_alpha4_arr*differences) - t8 = np.cumsum(i_alpha6_arr*differences) - h_arr = np.logspace(np.log10(1./n), np.log10(1.0), hsteps) - max_i_vector = (np.floor((n*h_arr))-2.).astype(int) - - gamma_pos = (35./(16*h_arr))*t1[max_i_vector]\ - - (105./(16*(h_arr**3)))*t2[max_i_vector]\ - + (105./(16*(h_arr**5)))*t3[max_i_vector]\ - - (35./(16*(h_arr**7)))*t4[max_i_vector] - - q1 = (35./(16*h_arr))*t5[max_i_vector]\ - + (105./(16*(h_arr**5)))*t7[max_i_vector]\ - - (105./(16*(h_arr**3)))*t6[max_i_vector]\ - - (35./(16*(h_arr**7)))*t8[max_i_vector] - - q2 = (35.*(1+alpha)/(16*h_arr))*t5[max_i_vector] \ - + (105.*(5+alpha)/(16*(h_arr**5)))*t7[max_i_vector] \ - - (105.*(3+alpha)/(16*(h_arr**3)))*t6[max_i_vector] \ - - (35.*(7+alpha)/(16*(h_arr**7)))*t8[max_i_vector] - - xi_arr = gamma_pos - 1. + q2/q1 + i_alpha2_arr = i_arr ** (2.0 + alpha) + i_alpha4_arr = i_arr ** (4.0 + alpha) + i_alpha6_arr = i_arr ** (6.0 + alpha) + t1 = np.cumsum(i_arr * differences) + t2 = np.cumsum(i3_arr * differences) + t3 = np.cumsum(i5_arr * differences) + t4 = np.cumsum(i7_arr * differences) + t5 = np.cumsum(i_alpha_arr * differences) + t6 = np.cumsum(i_alpha2_arr * differences) + t7 = np.cumsum(i_alpha4_arr * differences) + t8 = np.cumsum(i_alpha6_arr * differences) + h_arr = np.logspace(np.log10(1.0 / n), np.log10(1.0), hsteps) + max_i_vector = (np.floor(n * h_arr) - 2.0).astype(int) + + gamma_pos = ( + (35.0 / (16 * h_arr)) * t1[max_i_vector] + - (105.0 / (16 * (h_arr**3))) * t2[max_i_vector] + + (105.0 / (16 * (h_arr**5))) * t3[max_i_vector] + - (35.0 / (16 * (h_arr**7))) * t4[max_i_vector] + ) + + q1 = ( + (35.0 / (16 * h_arr)) * t5[max_i_vector] + + (105.0 / (16 * (h_arr**5))) * t7[max_i_vector] + - (105.0 / (16 * (h_arr**3))) * t6[max_i_vector] + - (35.0 / (16 * (h_arr**7))) * t8[max_i_vector] + ) + + q2 = ( + (35.0 * (1 + alpha) / (16 * h_arr)) * t5[max_i_vector] + + (105.0 * (5 + alpha) / (16 * (h_arr**5))) * t7[max_i_vector] + - (105.0 * (3 + alpha) / (16 * (h_arr**3))) * t6[max_i_vector] + - (35.0 * (7 + alpha) / (16 * (h_arr**7))) * t8[max_i_vector] + ) + + xi_arr = gamma_pos - 1.0 + q2 / q1 return h_arr, xi_arr -def kernel_type_dbs(ordered_data, hsteps, t_bootstrap = 0.5, - r_bootstrap = 500, alpha = 0.6, eps_stop = 1.0, - verbose = False, diagn_plots = False, base_seed=None): + +def kernel_type_dbs( + ordered_data, + hsteps, + t_bootstrap=0.5, + r_bootstrap=500, + alpha=0.6, + eps_stop=1.0, + verbose=False, + diagn_plots=False, + base_seed=None, +): """ - Function to perform double-bootstrap procedure for + Function to perform double-bootstrap procedure for moments estimator. Args: @@ -851,11 +1007,11 @@ def kernel_type_dbs(ordered_data, hsteps, t_bootstrap = 0.5, eps_stop: parameter controlling range of AMSE minimization. Defined as the fraction of order statistics to consider during the AMSE minimization step. - verbose: flag controlling bootstrap verbosity. + verbose: flag controlling bootstrap verbosity. diagn_plots: flag to switch on/off generation of AMSE diagnostic plots. base_seed: base random seed for reproducibility of bootstrap (default is None). - + Returns: h_star: fraction of order statistics optimal for estimation @@ -882,65 +1038,75 @@ def kernel_type_dbs(ordered_data, hsteps, t_bootstrap = 0.5, if verbose: logging.debug("Performing kernel double-bootstrap...") n = len(ordered_data) - eps_bootstrap = 0.5*(1+np.log(int(t_bootstrap*n))/np.log(n)) + eps_bootstrap = 0.5 * (1 + np.log(int(t_bootstrap * n)) / np.log(n)) + + base_rng = np.random.default_rng( + seed=base_seed + ) # Accept random seed for reproducibility. Default seed is None. - base_rng = np.random.default_rng(seed=base_seed) # Accept random seed for reproducibility. Default seed is None. - # first bootstrap with n1 sample size n1 = int(n**eps_bootstrap) samples_n1 = np.zeros(hsteps) good_counts1 = np.zeros(hsteps) - for i in range(r_bootstrap): + for _i in range(r_bootstrap): bs1_seed = base_rng.integers(0, 1_000_000) bs1_rng = np.random.default_rng(bs1_seed) - sample = bs1_rng.choice(ordered_data, n1, replace = True) + sample = bs1_rng.choice(ordered_data, n1, replace=True) sample[::-1].sort() _, xi2_arr = get_biweight_kernel_estimates(sample, hsteps, alpha) _, xi3_arr = get_triweight_kernel_estimates(sample, hsteps, alpha) - samples_n1 += (xi2_arr - xi3_arr)**2 - good_counts1[np.where((xi2_arr - xi3_arr)**2 != np.nan)] += 1 - max_index1 = (np.abs(np.logspace(np.log10(1./n1), np.log10(1.0), hsteps) - eps_stop)).argmin() - x1_arr = np.logspace(np.log10(1./n1), np.log10(1.0), hsteps) + samples_n1 += (xi2_arr - xi3_arr) ** 2 + good_counts1[np.where((xi2_arr - xi3_arr) ** 2 != np.nan)] += 1 + max_index1 = ( + np.abs(np.logspace(np.log10(1.0 / n1), np.log10(1.0), hsteps) - eps_stop) + ).argmin() + x1_arr = np.logspace(np.log10(1.0 / n1), np.log10(1.0), hsteps) averaged_delta = samples_n1 / good_counts1 - h1 = x1_arr[np.nanargmin(averaged_delta[:max_index1])] + h1 = x1_arr[np.nanargmin(averaged_delta[:max_index1])] if diagn_plots: n1_amse = averaged_delta - - + # second bootstrap with n2 sample size - n2 = int(n1*n1/float(n)) + n2 = int(n1 * n1 / float(n)) if n2 < hsteps: - sys.exit("Number of h points is larger than number "+\ - "of order statistics! Please either increase "+\ - "the size of 2nd bootstrap or decrease number "+\ - "of h grid points.") + sys.exit( + "Number of h points is larger than number " + + "of order statistics! Please either increase " + + "the size of 2nd bootstrap or decrease number " + + "of h grid points." + ) samples_n2 = np.zeros(hsteps) good_counts2 = np.zeros(hsteps) - for i in range(r_bootstrap): + for _i in range(r_bootstrap): bs2_seed = base_rng.integers(0, 1_000_000) bs2_rng = np.random.default_rng(bs2_seed) - sample = bs2_rng.choice(ordered_data, n2, replace = True) + sample = bs2_rng.choice(ordered_data, n2, replace=True) sample[::-1].sort() _, xi2_arr = get_biweight_kernel_estimates(sample, hsteps, alpha) _, xi3_arr = get_triweight_kernel_estimates(sample, hsteps, alpha) - samples_n2 += (xi2_arr - xi3_arr)**2 - good_counts2[np.where((xi2_arr - xi3_arr)**2 != np.nan)] += 1 - max_index2 = (np.abs(np.logspace(np.log10(1./n2), np.log10(1.0), hsteps) - eps_stop)).argmin() - x2_arr = np.logspace(np.log10(1./n2), np.log10(1.0), hsteps) + samples_n2 += (xi2_arr - xi3_arr) ** 2 + good_counts2[np.where((xi2_arr - xi3_arr) ** 2 != np.nan)] += 1 + max_index2 = ( + np.abs(np.logspace(np.log10(1.0 / n2), np.log10(1.0), hsteps) - eps_stop) + ).argmin() + x2_arr = np.logspace(np.log10(1.0 / n2), np.log10(1.0), hsteps) averaged_delta = samples_n2 / good_counts2 h2 = x2_arr[np.nanargmin(averaged_delta[:max_index2])] if diagn_plots: n2_amse = averaged_delta - - A = (143.*((np.log(n1) + np.log(h1))**2)/(3*(np.log(n1) - 13. * np.log(h1))**2))\ - **(-np.log(h1)/np.log(n1)) - - h_star = (h1*h1/float(h2)) * A + + A = ( + 143.0 + * ((np.log(n1) + np.log(h1)) ** 2) + / (3 * (np.log(n1) - 13.0 * np.log(h1)) ** 2) + ) ** (-np.log(h1) / np.log(n1)) + + h_star = (h1 * h1 / float(h2)) * A if h_star > 1: logging.warning("WARNING: estimated threshold is larger than the size of data!") logging.warning("WARNING: optimal h is set to 1...") - h_star = 1. + h_star = 1.0 if verbose: logging.info("--- Kernel-type double-bootstrap information ---") @@ -963,10 +1129,19 @@ def kernel_type_dbs(ordered_data, hsteps, t_bootstrap = 0.5, max_k_index2 = None return h_star, x1_arr, n1_amse, h1, max_k_index1, x2_arr, n2_amse, h2, max_k_index2 -def kernel_type_estimator(ordered_data, hsteps, alpha = 0.6, - bootstrap = True, t_bootstrap = 0.5, - r_bootstrap = 500, verbose = False, - diagn_plots = False, eps_stop = 0.99, base_seed = None): + +def kernel_type_estimator( + ordered_data, + hsteps, + alpha=0.6, + bootstrap=True, + t_bootstrap=0.5, + r_bootstrap=500, + verbose=False, + diagn_plots=False, + eps_stop=0.99, + base_seed=None, +): """ Function to calculate kernel-type estimator for a given dataset. If bootstrap flag is True, double-bootstrap procedure @@ -989,7 +1164,7 @@ def kernel_type_estimator(ordered_data, hsteps, alpha = 0.6, eps_stop: parameter controlling range of AMSE minimization. Defined as the fraction of order statistics to consider during the AMSE minimization step. - verbose: flag controlling bootstrap verbosity. + verbose: flag controlling bootstrap verbosity. diagn_plots: flag to switch on/off generation of AMSE diagnostic plots. base_seed: base random seed for reproducibility of bootstrap (default is None). @@ -1010,47 +1185,76 @@ def kernel_type_estimator(ordered_data, hsteps, alpha = 0.6, """ n = len(ordered_data) - h_arr, xi_arr = get_biweight_kernel_estimates(ordered_data, hsteps, - alpha = alpha) + h_arr, xi_arr = get_biweight_kernel_estimates(ordered_data, hsteps, alpha=alpha) if bootstrap: - results = kernel_type_dbs(ordered_data, hsteps, - t_bootstrap = t_bootstrap, - alpha = alpha, r_bootstrap = r_bootstrap, - verbose = verbose, diagn_plots = diagn_plots, - eps_stop = eps_stop, base_seed = base_seed) - h_star, x1_arr, n1_amse, h1, max_index1, x2_arr, n2_amse, h2, max_index2 = results - while h_star == None: + results = kernel_type_dbs( + ordered_data, + hsteps, + t_bootstrap=t_bootstrap, + alpha=alpha, + r_bootstrap=r_bootstrap, + verbose=verbose, + diagn_plots=diagn_plots, + eps_stop=eps_stop, + base_seed=base_seed, + ) + h_star, x1_arr, n1_amse, h1, max_index1, x2_arr, n2_amse, h2, max_index2 = ( + results + ) + while h_star is None: logging.debug("Resampling...") - results = kernel_type_dbs(ordered_data, hsteps, - t_bootstrap = t_bootstrap, - alpha = alpha, r_bootstrap = r_bootstrap, - verbose = verbose, diagn_plots = diagn_plots, - eps_stop = eps_stop, base_seed = base_seed) - h_star, x1_arr, n1_amse, h1, max_index1, x2_arr, n2_amse, h2, max_index2 = results - - #get k index which corresponds to h_star + results = kernel_type_dbs( + ordered_data, + hsteps, + t_bootstrap=t_bootstrap, + alpha=alpha, + r_bootstrap=r_bootstrap, + verbose=verbose, + diagn_plots=diagn_plots, + eps_stop=eps_stop, + base_seed=base_seed, + ) + h_star, x1_arr, n1_amse, h1, max_index1, x2_arr, n2_amse, h2, max_index2 = ( + results + ) + + # get k index which corresponds to h_star k_star = np.argmin(np.abs(h_arr - h_star)) xi_star = xi_arr[k_star] k_arr = [] - k_star = int(np.floor(h_arr[k_star]*n))-1 - k_arr = np.floor((h_arr * n)) + k_star = int(np.floor(h_arr[k_star] * n)) - 1 + k_arr = np.floor(h_arr * n) if xi_star <= 0: logging.info("Kernel-type estimated gamma: infinity (xi <= 0).") else: - logging.info(f"Kernel-type estimated gamma: {1 + 1./xi_star}") + logging.info(f"Kernel-type estimated gamma: {1 + 1.0 / xi_star}") else: k_star, xi_star = None, None - x1_arr, n1_amse, h1, max_index1 = 4*[None] - x2_arr, n2_amse, h2, max_index2 = 4*[None] + x1_arr, n1_amse, h1, max_index1 = 4 * [None] + x2_arr, n2_amse, h2, max_index2 = 4 * [None] k_arr = np.floor(h_arr * n) - results = [np.array(k_arr), xi_arr, k_star, xi_star, x1_arr, n1_amse, h1, max_index1,\ - x2_arr, n2_amse, h2, max_index2] + results = [ + np.array(k_arr), + xi_arr, + k_star, + xi_star, + x1_arr, + n1_amse, + h1, + max_index1, + x2_arr, + n2_amse, + h2, + max_index2, + ] return results + # ==================================================== # ========== Pickands Tail Index Estimation ========== # ==================================================== + def pickands_estimator(ordered_data): """ Function to calculate Pickands estimator for the tail index. @@ -1067,12 +1271,12 @@ def pickands_estimator(ordered_data): to k-order statistics provided in k_arr. """ n = len(ordered_data) - indices_k = np.arange(1, int(np.floor(n/4.))+1) - indices_2k = 2*indices_k - indices_4k = 4*indices_k - Z_k = ordered_data[indices_k-1] - Z_2k = ordered_data[indices_2k-1] - Z_4k = ordered_data[indices_4k-1] - xi_arr = (1./np.log(2)) * np.log((Z_k - Z_2k) / (Z_2k - Z_4k)) - k_arr = np.array([float(i) for i in range(1, int(np.floor(n/4.))+1)]) - return k_arr, xi_arr \ No newline at end of file + indices_k = np.arange(1, int(np.floor(n / 4.0)) + 1) + indices_2k = 2 * indices_k + indices_4k = 4 * indices_k + Z_k = ordered_data[indices_k - 1] + Z_2k = ordered_data[indices_2k - 1] + Z_4k = ordered_data[indices_4k - 1] + xi_arr = (1.0 / np.log(2)) * np.log((Z_k - Z_2k) / (Z_2k - Z_4k)) + k_arr = np.array([float(i) for i in range(1, int(np.floor(n / 4.0)) + 1)]) + return k_arr, xi_arr diff --git a/tests/test_datasets.py b/tests/test_datasets.py index 53bdcdd..95786a7 100644 --- a/tests/test_datasets.py +++ b/tests/test_datasets.py @@ -1,104 +1,116 @@ import os -import pytest -import numpy as np import tempfile + +import numpy as np +import pytest + from tailestim.datasets import TailData + def test_load_existing_data(): """Test loading an existing dataset""" # Test with CAIDA dataset data = TailData(name="CAIDA_KONECT") assert isinstance(data.data, np.ndarray) assert len(data.data) > 0 - + # Test with Libimseti dataset data = TailData(name="Libimseti_in_KONECT") assert isinstance(data.data, np.ndarray) assert len(data.data) > 0 - + # Test with Pareto dataset data = TailData(name="Pareto") assert isinstance(data.data, np.ndarray) assert len(data.data) > 0 + def test_nonexistent_data(): """Test handling of non-existent dataset""" with pytest.raises(FileNotFoundError): TailData(name="nonexistent_dataset") - + + def test_missing_parameters(): """Test handling of missing parameters""" with pytest.raises(ValueError): TailData() + def test_data_format(): """Test if data is properly formatted""" data = TailData("CAIDA_KONECT") - + # Data should be a numpy array assert isinstance(data.data, np.ndarray) - + # All values should be numeric assert np.issubdtype(data.data.dtype, np.number) - + + def test_representation(): """Test string representation of TailData""" data = TailData("CAIDA_KONECT") repr_str = repr(data) - + # Check if representation contains the name assert "CAIDA_KONECT" in repr_str - + # Check if representation contains the data length assert str(len(data.data)) in repr_str - + # Check format assert repr_str.startswith("TailData(") assert repr_str.endswith(")") + + def test_data_consistency(): """Test if loaded data is consistent with file content""" data = TailData(name="CAIDA_KONECT") - + # Get the data file path - examples_dir = os.path.join(os.path.dirname(os.path.dirname(__file__)), 'src', 'tailestim', 'data') - file_path = os.path.join(examples_dir, 'CAIDA_KONECT.dat') - + examples_dir = os.path.join( + os.path.dirname(os.path.dirname(__file__)), "src", "tailestim", "data" + ) + file_path = os.path.join(examples_dir, "CAIDA_KONECT.dat") + # Read the file manually to verify data - with open(file_path, 'r') as file: + with open(file_path) as file: lines = file.readlines() - + # Calculate total count from file total_count = sum(int(line.strip().split()[1]) for line in lines) - + # Verify the loaded data length matches the file content assert len(data.data) == total_count - + # Verify some values from the file match the loaded data first_line = lines[0].strip().split() first_value = float(first_line[0]) first_count = int(first_line[1]) assert np.all(data.data[:first_count] == first_value) - + + def test_data_load_custom_path(): """Test loading data from a custom path""" # Create a temporary file with test data - with tempfile.NamedTemporaryFile(mode='w', delete=False) as temp_file: + with tempfile.NamedTemporaryFile(mode="w", delete=False) as temp_file: temp_file.write("10 3\n") temp_file.write("20 2\n") temp_file.write("30 1\n") temp_path = temp_file.name - + try: # Load data from the custom path data = TailData(path=temp_path) - + # Verify the data was loaded correctly assert isinstance(data.data, np.ndarray) assert len(data.data) == 6 # 3 + 2 + 1 = 6 assert np.all(data.data[:3] == 10) assert np.all(data.data[3:5] == 20) assert data.data[5] == 30 - + # Check the string representation repr_str = repr(data) assert temp_path in repr_str @@ -106,21 +118,22 @@ def test_data_load_custom_path(): finally: # Clean up the temporary file os.unlink(temp_path) - + + def test_path_precedence(): """Test that path takes precedence over name when both are provided""" # Create a temporary file with test data - with tempfile.NamedTemporaryFile(mode='w', delete=False) as temp_file: + with tempfile.NamedTemporaryFile(mode="w", delete=False) as temp_file: temp_file.write("100 1\n") temp_path = temp_file.name - + try: # Load data with both name and path data = TailData(name="CAIDA_KONECT", path=temp_path) - + # Verify that the path was used, not the name assert len(data.data) == 1 assert data.data[0] == 100 finally: # Clean up the temporary file - os.unlink(temp_path) \ No newline at end of file + os.unlink(temp_path) diff --git a/tests/test_estimator_set.py b/tests/test_estimator_set.py index 8960a50..16ed425 100644 --- a/tests/test_estimator_set.py +++ b/tests/test_estimator_set.py @@ -1,12 +1,19 @@ +import matplotlib.pyplot as plt import numpy as np import pytest -import matplotlib.pyplot as plt + +from tailestim.datasets import TailData +from tailestim.estimators.estimator_set import TailEstimatorSet + pytestmark = [ - pytest.mark.filterwarnings("ignore:invalid value encountered in divide:RuntimeWarning"), - pytest.mark.filterwarnings("ignore:divide by zero encountered in divide:RuntimeWarning") + pytest.mark.filterwarnings( + "ignore:invalid value encountered in divide:RuntimeWarning" + ), + pytest.mark.filterwarnings( + "ignore:divide by zero encountered in divide:RuntimeWarning" + ), ] -from tailestim.estimators.estimator_set import TailEstimatorSet -from tailestim.datasets import TailData + def test_tail_estimator_set_initialization(): """Test that TailEstimatorSet can be initialized without data.""" @@ -17,58 +24,63 @@ def test_tail_estimator_set_initialization(): assert estimator_set.fig is None assert estimator_set.axes is None + def test_tail_estimator_set_fit(): """Test that TailEstimatorSet can fit data.""" # Generate Pareto distributed data np.random.seed(42) size = 1000 data = np.random.pareto(2, size) - + # Initialize estimator_set = TailEstimatorSet() - + # Fit data estimator_set.fit(data) - + # Check that data was stored and processed assert estimator_set.data is not None assert estimator_set.ordered_data is not None assert estimator_set.results is not None assert len(estimator_set.data) == size - assert estimator_set.ordered_data[0] >= estimator_set.ordered_data[-1] # Check ordering + assert ( + estimator_set.ordered_data[0] >= estimator_set.ordered_data[-1] + ) # Check ordering + def test_tail_estimator_set_plot(): """Test that TailEstimatorSet can generate plots.""" # Generate Pareto distributed data np.random.seed(42) data = np.random.pareto(2, 1000) - + estimator_set = TailEstimatorSet() estimator_set.fit(data) - + # Generate plots fig, axes = estimator_set.plot() - + # Check that plots were generated assert fig is not None assert axes is not None assert isinstance(fig, plt.Figure) assert isinstance(axes, np.ndarray) assert axes.shape == (3, 2) # 3 rows, 2 columns of plots - + # Check that the figure and axes are stored in the object assert estimator_set.fig is fig assert estimator_set.axes is axes - + # Clean up plt.close(fig) + def test_tail_estimator_set_diagnostic_plot(): """Test that TailEstimatorSet can generate diagnostic plots.""" # Generate Pareto distributed data np.random.seed(42) data = np.random.pareto(2, 1000) - + # Initialize with data and enable diagnostic plots estimator_set = TailEstimatorSet( bootstrap_flag=True, @@ -77,103 +89,106 @@ def test_tail_estimator_set_diagnostic_plot(): ) estimator_set.fit(data) - + # Generate diagnostic plots fig_d, axes_d = estimator_set.plot_diagnostics() - + # Check that diagnostic plots were generated assert fig_d is not None assert axes_d is not None assert isinstance(fig_d, plt.Figure) assert isinstance(axes_d, np.ndarray) - + # Clean up plt.close(fig_d) + def test_tail_estimator_set_with_built_in_dataset(): """Test TailEstimatorSet with a built-in dataset.""" try: # Load a built-in dataset - data = TailData(name='CAIDA_KONECT').data - + data = TailData(name="CAIDA_KONECT").data + estimator_set = TailEstimatorSet( bootstrap_flag=True, diagnostic_plots=True, - r_bootstrap=100 # Reduce bootstrap iterations for faster testing + r_bootstrap=100, # Reduce bootstrap iterations for faster testing ) estimator_set.fit(data) - + # Generate plots fig, axes = estimator_set.plot() assert fig is not None assert axes is not None plt.close(fig) - + # Generate diagnostic plots fig_d, axes_d = estimator_set.plot_diagnostics() assert fig_d is not None assert axes_d is not None plt.close(fig_d) - + except FileNotFoundError: pytest.skip("Built-in dataset not found, skipping test") + def test_tail_estimator_set_errors(): """Test that TailEstimatorSet raises appropriate errors.""" # Initialize without data estimator_set = TailEstimatorSet() - + # Attempt to plot without fitting data with pytest.raises(ValueError, match="No data has been fitted"): estimator_set.plot() - + # Attempt to get diagnostic plots without fitting data with pytest.raises(ValueError, match="No data has been fitted"): estimator_set.plot_diagnostics() - + # Fit data but with bootstrap disabled np.random.seed(42) data = np.random.pareto(2, 1000) estimator_set = TailEstimatorSet(bootstrap_flag=False) estimator_set.fit(data) - + # Attempt to get diagnostic plots with bootstrap disabled - with pytest.raises(ValueError, match="Diagnostic plots require bootstrap to be enabled"): + with pytest.raises( + ValueError, match="Diagnostic plots require bootstrap to be enabled" + ): estimator_set.plot_diagnostics() - + # Fit data with bootstrap enabled but diagnostic plots disabled estimator_set = TailEstimatorSet(bootstrap_flag=True, diagnostic_plots=False) estimator_set.fit(data) - + # Attempt to get diagnostic plots with diagnostic plots disabled with pytest.raises(ValueError, match="Diagnostic plots are not enabled"): estimator_set.plot_diagnostics() + def test_tail_estimator_set_parameters(): """Test that TailEstimatorSet parameters are correctly stored and retrieved.""" # Generate data np.random.seed(42) size = 1000 data = np.random.pareto(2, size) - + # Initialize with custom parameters custom_bins = 50 custom_r_smooth = 3 custom_alpha = 0.7 estimator_set = TailEstimatorSet( - number_of_bins=custom_bins, - r_smooth=custom_r_smooth, - alpha=custom_alpha + number_of_bins=custom_bins, r_smooth=custom_r_smooth, alpha=custom_alpha ) estimator_set.fit(data) - + # Get parameters params = estimator_set.get_params() - + # Check that parameters were correctly stored - assert params['number_of_bins'] == custom_bins - assert params['r_smooth'] == custom_r_smooth - assert params['alpha'] == custom_alpha - assert params['data_length'] == size \ No newline at end of file + assert params["number_of_bins"] == custom_bins + assert params["r_smooth"] == custom_r_smooth + assert params["alpha"] == custom_alpha + assert params["data_length"] == size diff --git a/tests/test_tail_methods.py b/tests/test_tail_methods.py index 2e0b836..3618110 100644 --- a/tests/test_tail_methods.py +++ b/tests/test_tail_methods.py @@ -1,26 +1,35 @@ import numpy as np import pytest -pytestmark = [ - pytest.mark.filterwarnings("ignore:invalid value encountered in divide:RuntimeWarning"), - pytest.mark.filterwarnings("ignore:divide by zero encountered in divide:RuntimeWarning") -] + +from tailestim.estimators.hill import HillEstimator +from tailestim.estimators.kernel import KernelTypeEstimator +from tailestim.estimators.moments import MomentsEstimator +from tailestim.estimators.pickands import PickandsEstimator +from tailestim.estimators.smooth_hill import SmoothHillEstimator from tailestim.estimators.tail_methods import ( add_uniform_noise, + get_ccdf, get_distribution, - get_ccdf ) -from tailestim.estimators.hill import HillEstimator -from tailestim.estimators.smooth_hill import SmoothHillEstimator -from tailestim.estimators.moments import MomentsEstimator -from tailestim.estimators.kernel import KernelTypeEstimator -from tailestim.estimators.pickands import PickandsEstimator + +pytestmark = [ + pytest.mark.filterwarnings( + "ignore:invalid value encountered in divide:RuntimeWarning" + ), + pytest.mark.filterwarnings( + "ignore:divide by zero encountered in divide:RuntimeWarning" + ), +] + # Test data preprocessing functions def test_add_uniform_noise(): # Test with valid input data = np.array([1, 2, 3, 4, 5]) noisy_data = add_uniform_noise(data, p=1) - assert len(noisy_data) <= len(data) # May be shorter due to negative value filtering + assert len(noisy_data) <= len( + data + ) # May be shorter due to negative value filtering assert np.all(noisy_data > 0) # All values should be positive # Test with invalid p @@ -41,20 +50,27 @@ def test_add_uniform_noise(): assert len(result4) > 0 # Should return valid array assert np.all(result4 > 0) # All values should be positive + def test_get_distribution(): data = np.array([1, 2, 2, 3, 3, 3, 4, 4, 5]) x, y = get_distribution(data, number_of_bins=5) assert len(x) == len(y) assert np.all(np.array(y) >= 0) # PDF values should be non-negative + def test_get_ccdf(): data = np.array([1, 2, 2, 3, 3, 3, 4, 4, 5]) # uniques are returned in descending order # ccdf is returned for each value in uniques (in descending order) uniques, ccdf = get_ccdf(data) assert len(uniques) == len(ccdf) - assert np.all(ccdf >= 0) and np.all(ccdf <= 1) # CCDF values should be between 0 and 1 - assert ccdf[0] == 0 # The first element of returned ccdf object is CCDF for last unique degree + assert np.all(ccdf >= 0) and np.all( + ccdf <= 1 + ) # CCDF values should be between 0 and 1 + assert ( + ccdf[0] == 0 + ) # The first element of returned ccdf object is CCDF for last unique degree + # Test Hill estimator def test_hill_estimator(): @@ -62,35 +78,35 @@ def test_hill_estimator(): np.random.seed(42) alpha = 2.0 size = 1000 - data = (1/np.random.uniform(0, 1, size))**(1/alpha) - + data = (1 / np.random.uniform(0, 1, size)) ** (1 / alpha) + # Test without bootstrap estimator = HillEstimator(bootstrap=False) estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'k_arr_') - assert hasattr(res, 'xi_arr_') + assert hasattr(res, "k_arr_") + assert hasattr(res, "xi_arr_") assert len(res.k_arr_) == len(res.xi_arr_) # Test that params are returned params = estimator.get_params() assert params is not None - + # Test with bootstrap estimator = HillEstimator(bootstrap=True, r_bootstrap=100) estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'estimator') + assert hasattr(res, "estimator") assert isinstance(res.estimator, HillEstimator) - assert hasattr(res, 'k_star_') - assert hasattr(res, 'xi_star_') - assert hasattr(res, 'gamma_') + assert hasattr(res, "k_star_") + assert hasattr(res, "xi_star_") + assert hasattr(res, "gamma_") assert res.gamma_ is not None # Test with bootstrap with seed. Run multiple times to ensure consistency. bs1_results = [] bs2_results = [] - for i in range(3): + for _i in range(3): estimator = HillEstimator(bootstrap=True, base_seed=42, r_bootstrap=100) estimator.fit(data) res = estimator.get_result() @@ -99,37 +115,155 @@ def test_hill_estimator(): bs1_results.append(bs1_kmin) bs2_results.append(bs2_kmin) # Assert that all values in bs1_results are the same - assert all(x == bs1_results[0] for x in bs1_results), "Bootstrap results with same seed should be identical" - assert all(x == bs2_results[0] for x in bs2_results), "Second bootstrap results with same seed should be identical" + assert all(x == bs1_results[0] for x in bs1_results), ( + "Bootstrap results with same seed should be identical" + ) + assert all(x == bs2_results[0] for x in bs2_results), ( + "Second bootstrap results with same seed should be identical" + ) + def test_hill_estimator_max_resample(): """Test that HillEstimator raises RuntimeError when max_resample is exceeded.""" # Near-constant degree sequence (n=115). # With n < 200, the boundary shift int(0.005*n) == 0, so k2 > k1 # persists indefinitely and max_resample is always hit. - data = np.array([ - 12,12,12,12,11,12,12,12,11,11,10,10,10,11,10,12,11,11,11,11, - 11,11,11,11,10,11,10,11,9,11,11,11,11,10,11,11,8,11,11,11, - 11,10,7,11,11,11,11,11,11,11,9,11,10,12,10,11,10,10,10,8, - 11,11,11,9,11,11,11,12,11,11,11,10,11,11,11,10,11,11,11,11, - 11,11,11,11,11,9,11,11,12,11,9,11,11,10,10,10,10,8,11,10, - 11,10,10,10,12,10,11,10,10,11,11,11,10,10,11, - ], dtype=float) + data = np.array( + [ + 12, + 12, + 12, + 12, + 11, + 12, + 12, + 12, + 11, + 11, + 10, + 10, + 10, + 11, + 10, + 12, + 11, + 11, + 11, + 11, + 11, + 11, + 11, + 11, + 10, + 11, + 10, + 11, + 9, + 11, + 11, + 11, + 11, + 10, + 11, + 11, + 8, + 11, + 11, + 11, + 11, + 10, + 7, + 11, + 11, + 11, + 11, + 11, + 11, + 11, + 9, + 11, + 10, + 12, + 10, + 11, + 10, + 10, + 10, + 8, + 11, + 11, + 11, + 9, + 11, + 11, + 11, + 12, + 11, + 11, + 11, + 10, + 11, + 11, + 11, + 10, + 11, + 11, + 11, + 11, + 11, + 11, + 11, + 11, + 11, + 9, + 11, + 11, + 12, + 11, + 9, + 11, + 11, + 10, + 10, + 10, + 10, + 8, + 11, + 10, + 11, + 10, + 10, + 10, + 12, + 10, + 11, + 10, + 10, + 11, + 11, + 11, + 10, + 10, + 11, + ], + dtype=float, + ) estimator = HillEstimator(bootstrap=True, max_resample=3, base_seed=42) with pytest.raises(RuntimeError, match="failed to converge after 3 resampling"): estimator.fit(data) + def test_smooth_hill_estimator(): np.random.seed(42) data = np.random.pareto(2, 1000) - + estimator = SmoothHillEstimator(r_smooth=2) estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'estimator') + assert hasattr(res, "estimator") assert isinstance(res.estimator, SmoothHillEstimator) - assert hasattr(res, 'k_arr_') - assert hasattr(res, 'xi_arr_') + assert hasattr(res, "k_arr_") + assert hasattr(res, "xi_arr_") assert len(res.k_arr_) == len(res.xi_arr_) assert np.all(np.isfinite(res.xi_arr_)) @@ -137,38 +271,39 @@ def test_smooth_hill_estimator(): params = estimator.get_params() assert params is not None + # Test moments estimator def test_moments_estimator(): np.random.seed(42) data = np.random.pareto(2, 1000) - + # Test without bootstrap estimator = MomentsEstimator(bootstrap=False) estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'k_arr_') - assert hasattr(res, 'xi_arr_') + assert hasattr(res, "k_arr_") + assert hasattr(res, "xi_arr_") assert len(res.k_arr_) == len(res.xi_arr_) # Test that params are returned params = estimator.get_params() assert params is not None - + # Test with bootstrap estimator = MomentsEstimator(bootstrap=True, r_bootstrap=100) estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'estimator') + assert hasattr(res, "estimator") assert isinstance(res.estimator, MomentsEstimator) - assert hasattr(res, 'k_star_') - assert hasattr(res, 'xi_star_') + assert hasattr(res, "k_star_") + assert hasattr(res, "xi_star_") assert res.k_star_ is not None assert res.xi_star_ is not None # Test with bootstrap with seed. Run multiple times to ensure consistency. bs1_results = [] bs2_results = [] - for i in range(3): + for _i in range(3): estimator = MomentsEstimator(bootstrap=True, base_seed=42, r_bootstrap=100) estimator.fit(data) res = estimator.get_result() @@ -177,42 +312,49 @@ def test_moments_estimator(): bs1_results.append(bs1_kmin) bs2_results.append(bs2_kmin) # Assert that all values in bs1_results are the same - assert all(x == bs1_results[0] for x in bs1_results), "Bootstrap results with same seed should be identical" - assert all(x == bs2_results[0] for x in bs2_results), "Second bootstrap results with same seed should be identical" + assert all(x == bs1_results[0] for x in bs1_results), ( + "Bootstrap results with same seed should be identical" + ) + assert all(x == bs2_results[0] for x in bs2_results), ( + "Second bootstrap results with same seed should be identical" + ) + # Test kernel estimator def test_kernel_type_estimator(): np.random.seed(42) data = np.random.pareto(2, 1000) - + # Test without bootstrap estimator = KernelTypeEstimator(hsteps=50, bootstrap=False) estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'k_arr_') - assert hasattr(res, 'xi_arr_') + assert hasattr(res, "k_arr_") + assert hasattr(res, "xi_arr_") assert len(res.k_arr_) == len(res.xi_arr_) # Test that params are returned params = estimator.get_params() assert params is not None - + # Test with bootstrap estimator = KernelTypeEstimator(hsteps=50, bootstrap=True, r_bootstrap=100) estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'estimator') + assert hasattr(res, "estimator") assert isinstance(res.estimator, KernelTypeEstimator) - assert hasattr(res, 'k_star_') - assert hasattr(res, 'xi_star_') + assert hasattr(res, "k_star_") + assert hasattr(res, "xi_star_") assert res.k_star_ is not None assert res.xi_star_ is not None # Test with bootstrap with seed. Run multiple times to ensure consistency. bs1_results = [] bs2_results = [] - for i in range(3): - estimator = KernelTypeEstimator(hsteps=50, bootstrap=True, base_seed=42, r_bootstrap=100) + for _i in range(3): + estimator = KernelTypeEstimator( + hsteps=50, bootstrap=True, base_seed=42, r_bootstrap=100 + ) estimator.fit(data) res = estimator.get_result() bs1_kmin = res.bootstrap_results_.first_bootstrap_.h_min_ @@ -220,24 +362,31 @@ def test_kernel_type_estimator(): bs1_results.append(bs1_kmin) bs2_results.append(bs2_kmin) # Assert that all values in bs1_results are the same - assert all(x == bs1_results[0] for x in bs1_results), "Bootstrap results with same seed should be identical" - assert all(x == bs2_results[0] for x in bs2_results), "Second bootstrap results with same seed should be identical" + assert all(x == bs1_results[0] for x in bs1_results), ( + "Bootstrap results with same seed should be identical" + ) + assert all(x == bs2_results[0] for x in bs2_results), ( + "Second bootstrap results with same seed should be identical" + ) + # Test Pickands estimator def test_pickands_estimator(): np.random.seed(42) data = np.random.pareto(2, 1000) - + estimator = PickandsEstimator() estimator.fit(data) res = estimator.get_result() - assert hasattr(res, 'estimator') + assert hasattr(res, "estimator") assert isinstance(res.estimator, PickandsEstimator) - assert hasattr(res, 'k_arr_') - assert hasattr(res, 'xi_arr_') + assert hasattr(res, "k_arr_") + assert hasattr(res, "xi_arr_") assert len(res.k_arr_) == len(res.xi_arr_) - assert len(res.k_arr_) <= len(data) // 4 # Pickands can only estimate up to n/4 order statistics + assert ( + len(res.k_arr_) <= len(data) // 4 + ) # Pickands can only estimate up to n/4 order statistics # Test that params are returned params = estimator.get_params() - assert params is not None \ No newline at end of file + assert params is not None diff --git a/tests/test_tailestimation_validation.py b/tests/test_tailestimation_validation.py index e98ecb7..3236a21 100644 --- a/tests/test_tailestimation_validation.py +++ b/tests/test_tailestimation_validation.py @@ -12,11 +12,12 @@ - Medium/Large: Synthetic Pareto distributions (1,000 and 5,000 points) """ +import importlib.util import sys +from pathlib import Path + import numpy as np import pytest -from pathlib import Path -import importlib.util # Add tail-estimation to path TAIL_EST_PATH = Path(__file__).parent.parent / "tail-estimation" / "Python3" @@ -667,7 +668,7 @@ def test_pdf_ccdf_data(self, dataset_name, dataset_func): x_ccdf_old, y_ccdf_old = tail_estimation_old.get_ccdf(ordered_data) # tailestim - from tailestim.estimators.tail_methods import get_distribution, get_ccdf + from tailestim.estimators.tail_methods import get_ccdf, get_distribution x_pdf_new, y_pdf_new = get_distribution(ordered_data, number_of_bins=30) x_ccdf_new, y_ccdf_new = get_ccdf(ordered_data) From b2d138747b59e8072769990d872f9e29dd847f1e Mon Sep 17 00:00:00 2001 From: Minami Ueda Date: Fri, 27 Mar 2026 11:44:45 -0400 Subject: [PATCH 3/3] Add RUF, NPY, C4 ruff rules and fix violations --- pyproject.toml | 8 +++++++- src/tailestim/__init__.py | 4 ++-- src/tailestim/estimators/__init__.py | 2 +- src/tailestim/estimators/result.py | 6 +++--- tests/test_tailestimation_validation.py | 10 ++++------ 5 files changed, 17 insertions(+), 13 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 59266ec..243cf85 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -94,13 +94,19 @@ select = [ "I", # isort "UP", # pyupgrade "B", # flake8-bugbear + "RUF", # ruff-specific rules + "NPY", # numpy-specific rules + "C4", # flake8-comprehensions ] ignore = [ - "E501", # line length handled by ruff format + "E501", # line length handled by ruff format + "RUF001", # ambiguous unicode in strings (intentional math symbols: γ, ξ) + "RUF002", # ambiguous unicode in docstrings ] [tool.ruff.lint.per-file-ignores] "src/tailestim/estimators/plot/*" = ["F841", "E741"] +"tests/*" = ["NPY002"] [tool.ruff.lint.isort] known-first-party = ["tailestim"] diff --git a/src/tailestim/__init__.py b/src/tailestim/__init__.py index 7c2c67c..20f0b0b 100644 --- a/src/tailestim/__init__.py +++ b/src/tailestim/__init__.py @@ -12,12 +12,12 @@ from .estimators.smooth_hill import SmoothHillEstimator __all__ = [ - "TailData", "BaseTailEstimator", "HillEstimator", - "MomentsEstimator", "KernelTypeEstimator", + "MomentsEstimator", "PickandsEstimator", "SmoothHillEstimator", + "TailData", "TailEstimatorSet", ] diff --git a/src/tailestim/estimators/__init__.py b/src/tailestim/estimators/__init__.py index a7023f9..851a59a 100644 --- a/src/tailestim/estimators/__init__.py +++ b/src/tailestim/estimators/__init__.py @@ -11,8 +11,8 @@ __all__ = [ "BaseTailEstimator", "HillEstimator", - "MomentsEstimator", "KernelTypeEstimator", + "MomentsEstimator", "PickandsEstimator", "SmoothHillEstimator", "TailEstimatorSet", diff --git a/src/tailestim/estimators/result.py b/src/tailestim/estimators/result.py index c034428..531b20c 100644 --- a/src/tailestim/estimators/result.py +++ b/src/tailestim/estimators/result.py @@ -1,4 +1,4 @@ -from typing import Any, Dict +from typing import Any, ClassVar, Dict, Optional import numpy as np @@ -34,7 +34,7 @@ class TailEstimatorResult: """ # Mapping of keys to human-readable labels - _key_labels = { + _key_labels: ClassVar[Dict[str, str]] = { "estimator": "Estimator", "xi_star_": "Tail index (ξ)", "gamma_": "Power law exponent (γ)", @@ -50,7 +50,7 @@ class TailEstimatorResult: "x_arr_": "Fraction of order statistics", } - def __init__(self, initial_data: Dict[str, Any] = None, **kwargs): + def __init__(self, initial_data: Optional[Dict[str, Any]] = None, **kwargs): data = initial_data or {} data.update(kwargs) # Recursively convert nested dicts to TailEstimatorResult objects. diff --git a/tests/test_tailestimation_validation.py b/tests/test_tailestimation_validation.py index 3236a21..6d50695 100644 --- a/tests/test_tailestimation_validation.py +++ b/tests/test_tailestimation_validation.py @@ -123,9 +123,7 @@ def generate_large_pareto(seed=456): ] # Full test suite including large datasets -FULL_DATASETS = DATASET_CONFIGS + [ - ("large", DatasetGenerator.generate_large_pareto), -] +FULL_DATASETS = [*DATASET_CONFIGS, ("large", DatasetGenerator.generate_large_pareto)] # ============================================================================ @@ -238,7 +236,7 @@ def test_hill_basic_no_bootstrap(self, dataset_name, dataset_func): # tail-estimation ordered_data = np.sort(data)[::-1] results_old = tail_estimation_old.hill_estimator(ordered_data, bootstrap=False) - k_arr_old, xi_arr_old, k_star_old, xi_star_old = results_old[:4] + k_arr_old, xi_arr_old, _k_star_old, _xi_star_old = results_old[:4] # tailestim from tailestim.estimators import HillEstimator @@ -479,7 +477,7 @@ class TestPickandsEstimator: def test_pickands_equivalence(self, dataset): """Test Pickands estimator (no bootstrap) across datasets.""" - data, dataset_name, description = dataset + data, _dataset_name, description = dataset # tail-estimation ordered_data = np.sort(data)[::-1] @@ -565,7 +563,7 @@ def test_all_estimators_on_dataset(self, dataset_name, dataset_func): kernel_succeeded = True except ValueError as e: # Kernel can fail on degenerate data (e.g., all values identical) - print(f"Kernel: Failed ({str(e)})") + print(f"Kernel: Failed ({e!s})") kernel_succeeded = False # Test Pickands