From 0b1b07e2f0b1b812f9c98d587554ae0146f13e08 Mon Sep 17 00:00:00 2001 From: Maani Beigy Date: Fri, 16 Aug 2024 11:37:35 +0200 Subject: [PATCH] :sparkles: adds conf_limits_nct_minimize_scalar --- .logs/bandit.json | 34 ++- .logs/bandit.txt | 4 +- .logs/black.txt | 2 +- .logs/complexity.txt | Bin 2159 -> 2285 bytes .logs/docstring.txt | Bin 3081 -> 3243 bytes .logs/maintainability.json | 2 +- .logs/maintainability.txt | Bin 817 -> 889 bytes .logs/mypy.txt | Bin 93 -> 93 bytes .logs/pylint-log.txt | Bin 287 -> 291 bytes CONTRIBUTING.md | 4 +- assets/docs/pycvcqv/examples.ipynb | 398 ++++++++++++++++++++++++++--- poetry.lock | 44 +++- pycvcqv/noncentralt.py | 112 ++++++++ pycvcqv/sanitizers.py | 31 ++- pyproject.toml | 1 + requirements.txt | Bin 845 -> 981 bytes tests/test_checkers.py | 17 +- tests/test_noncentralt.py | 53 ++++ 18 files changed, 630 insertions(+), 72 deletions(-) create mode 100644 pycvcqv/noncentralt.py create mode 100644 tests/test_noncentralt.py diff --git a/.logs/bandit.json b/.logs/bandit.json index f1d8dd0..a6451a8 100644 --- a/.logs/bandit.json +++ b/.logs/bandit.json @@ -1,6 +1,6 @@ { "errors": [], - "generated_at": "2024-08-16T06:33:16Z", + "generated_at": "2024-08-16T09:31:50Z", "metrics": { "_totals": { "CONFIDENCE.HIGH": 0, @@ -11,7 +11,7 @@ "SEVERITY.LOW": 0, "SEVERITY.MEDIUM": 0, "SEVERITY.UNDEFINED": 0, - "loc": 1609, + "loc": 1740, "nosec": 0, "skipped_tests": 0 }, @@ -106,6 +106,19 @@ "nosec": 0, "skipped_tests": 0 }, + "pycvcqv\\noncentralt.py": { + "CONFIDENCE.HIGH": 0, + "CONFIDENCE.LOW": 0, + "CONFIDENCE.MEDIUM": 0, + "CONFIDENCE.UNDEFINED": 0, + "SEVERITY.HIGH": 0, + "SEVERITY.LOW": 0, + "SEVERITY.MEDIUM": 0, + "SEVERITY.UNDEFINED": 0, + "loc": 94, + "nosec": 0, + "skipped_tests": 0 + }, "pycvcqv\\prepare_output.py": { "CONFIDENCE.HIGH": 0, "CONFIDENCE.LOW": 0, @@ -128,7 +141,7 @@ "SEVERITY.LOW": 0, "SEVERITY.MEDIUM": 0, "SEVERITY.UNDEFINED": 0, - "loc": 31, + "loc": 35, "nosec": 0, "skipped_tests": 0 }, @@ -193,7 +206,7 @@ "SEVERITY.LOW": 0, "SEVERITY.MEDIUM": 0, "SEVERITY.UNDEFINED": 0, - "loc": 123, + "loc": 120, "nosec": 0, "skipped_tests": 0 }, @@ -222,6 +235,19 @@ "loc": 231, "nosec": 0, "skipped_tests": 0 + }, + "tests\\test_noncentralt.py": { + "CONFIDENCE.HIGH": 0, + "CONFIDENCE.LOW": 0, + "CONFIDENCE.MEDIUM": 0, + "CONFIDENCE.UNDEFINED": 0, + "SEVERITY.HIGH": 0, + "SEVERITY.LOW": 0, + "SEVERITY.MEDIUM": 0, + "SEVERITY.UNDEFINED": 0, + "loc": 36, + "nosec": 0, + "skipped_tests": 0 } }, "results": [] diff --git a/.logs/bandit.txt b/.logs/bandit.txt index 34bda66..9e61d34 100644 --- a/.logs/bandit.txt +++ b/.logs/bandit.txt @@ -1,10 +1,10 @@ -Run started:2024-08-16 06:33:18.488005 +Run started:2024-08-16 09:31:52.510266 Test results: No issues identified. Code scanned: - Total lines of code: 1609 + Total lines of code: 1740 Total lines skipped (#nosec): 0 Total potential issues skipped due to specifically being disabled (e.g., #nosec BXXX): 0 diff --git a/.logs/black.txt b/.logs/black.txt index cba5dd2..5f53f85 100644 --- a/.logs/black.txt +++ b/.logs/black.txt @@ -1,2 +1,2 @@ All done! f370 -17 files would be left unchanged. +19 files would be left unchanged. diff --git a/.logs/complexity.txt b/.logs/complexity.txt index fcf138a5e8241e1ea3e46bbc43ec6b6b89dbca80..64c3984eac600dadf3e358d454b851da457a680e 100644 GIT binary patch delta 149 zcmaDa@K$gG4~u9XLq0d(OV8vj-pa7HtDM(|82TEiz zfD>hCtGgO>^@f UHX#;9OQ3*0voV9^WJm530FoyXlK=n! delta 32 ocmZ22*(tH%58LK>>_RMzrjs{v8?u@*C@|Da-pDLG`4{&_0J6pkwEzGB diff --git a/.logs/maintainability.json b/.logs/maintainability.json index 49f4907..60b9140 100644 --- a/.logs/maintainability.json +++ b/.logs/maintainability.json @@ -1,3 +1,3 @@ { - "maintainability": "91.4%" + "maintainability": "90.1%" } diff --git a/.logs/maintainability.txt b/.logs/maintainability.txt index 185aeca402a9b073eb2e31a882497d76a7dcfd21..2ac3d1e06d8e9476216a571d6d701947401ea9f1 100644 GIT binary patch delta 58 zcmdnU_LFVHG)B=phJ1!RhGd3RAX&mt#E{64!%#9&QIy@7L65oBo00sz904oCn1 delta 27 jcmey#wvlbaG{(vI7%SK<81xtn8H^{c7vEgM#L5T&hRX;) diff --git a/.logs/mypy.txt b/.logs/mypy.txt index 60020b444e0763e29c16c77a3a5fede101d582ef..0f2069b3afffccc4527d211d19b8928e550209d4 100644 GIT binary patch delta 9 Qcma!zonXUgG0|2501m(ct^fc4 delta 9 Qcma!zonXUgHqllA01mhUtN;K2 diff --git a/.logs/pylint-log.txt b/.logs/pylint-log.txt index 909f57db3702ef93aa3b7c0b41b351c850adb365..f57ef2d478650a23efb5e46af8026a7df19de13f 100644 GIT binary patch delta 51 scmbQww3ul^&%~}gO+y9)20b7)V9;j(v2_>}7_{MHnhd-QTwu%v017+@DgXcg delta 49 ycmZ3?G@ofg&&0kQbxQ_421^E027LxY1_K5i1_cIfAgu?)rVN@4ybN4m%mo1bxCa#g diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index bc7768f..6d099d8 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -84,14 +84,14 @@ poetry run mypy --install-types --non-interactive pycvcqv/ tests/ 7. Unit tests and coverage ```powershell -poetry run pytest --cov-report term --cov pycvcqv tests/ +poetry run pytest --cov-report html --cov pycvcqv tests/ poetry run coverage-badge -o assets/images/coverage.svg -f ``` 8. Lint ```powershell -poetry run pylint pycvcqv +poetry run pylint pycvcqv tests ``` 9. Code-style check diff --git a/assets/docs/pycvcqv/examples.ipynb b/assets/docs/pycvcqv/examples.ipynb index 05b76a3..48790fb 100644 --- a/assets/docs/pycvcqv/examples.ipynb +++ b/assets/docs/pycvcqv/examples.ipynb @@ -11,7 +11,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## CV" + "## CV\n" ] }, { @@ -22,7 +22,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 1, @@ -31,33 +31,352 @@ } ], "source": [ - "from matplotlib import rcParams, cycler\n", + "from typing import Optional\n", + "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", + "from matplotlib import cycler, rcParams\n", + "from scipy.optimize import minimize, minimize_scalar\n", + "from scipy.stats import nct\n", "\n", "plt.ion()" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 113, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + " message: Solution found.\n", + " success: True\n", + " status: 0\n", + " fun: 2.28705961145802e-19\n", + " x: 4.815359140504376\n", + " nit: 13\n", + " nfev: 13" + ] + }, + "execution_count": 113, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ncp1 = 2.83\n", + "dof1 = 126\n", + "conf_level1 = 0.95\n", + "alpha_lower1 = (1 - conf_level1) / 2\n", + "min_ncp = min(-150, -5 * ncp1)\n", + "max_ncp = max(150, 5 * ncp1)\n", + "tol1 = 1e-9\n", + "\n", + "\n", + "def ci_nct_lower(val_of_interest):\n", + " return (nct.ppf(alpha_lower1, dof1, val_of_interest, loc=0) - ncp1) ** 2\n", + "\n", + "\n", + "minimize_scalar(\n", + " ci_nct_lower,\n", + " bounds=(min_ncp, max_ncp),\n", + " method=\"bounded\",\n", + " options={\"xatol\": tol1},\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 134, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4.815359140504376\n", + "0.8337502600175457\n" + ] + } + ], + "source": [ + "ncp1 = 2.83\n", + "dof1 = 126\n", + "conf_level1 = 0.95\n", + "alpha_lower1 = (1 - conf_level1) / 2\n", + "alpha_upper1 = (1 - conf_level1) / 2\n", + "min_ncp = min(-150, -5 * ncp1)\n", + "max_ncp = max(150, 5 * ncp1)\n", + "tol1 = 1e-9\n", + "\n", + "\n", + "def ci_nct_lower(val_of_interest):\n", + " return (nct.ppf(alpha_lower1, dof1, val_of_interest, loc=0) - ncp1) ** 2\n", + "\n", + "\n", + "def ci_nct_upper(val_of_interest):\n", + " return (nct.ppf(1 - alpha_upper1, dof1, val_of_interest, loc=0) - ncp1) ** 2\n", + "\n", + "\n", + "lower_limit = minimize_scalar(\n", + " ci_nct_lower,\n", + " bounds=(min_ncp, max_ncp),\n", + " method=\"bounded\",\n", + " options={\"xatol\": tol1},\n", + ")\n", + "\n", + "upper_limit = minimize_scalar(\n", + " ci_nct_upper,\n", + " bounds=(min_ncp, max_ncp),\n", + " method=\"bounded\",\n", + " options={\"xatol\": tol1},\n", + ")\n", + "print(lower_limit.x)\n", + "print(upper_limit.x)" + ] + }, + { + "cell_type": "code", + "execution_count": 135, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.024999999971943743" + ] + }, + "execution_count": 135, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nct.cdf(ncp1, dof1, lower_limit.x, loc=0)" + ] + }, + { + "cell_type": "code", + "execution_count": 138, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.02499999995262825" + ] + }, + "execution_count": 138, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "1 - nct.cdf(ncp1, dof1, upper_limit.x, loc=0)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", - "from scipy.optimize import minimize_scalar, minimize\n", - "from scipy.stats import nct\n", + "def conf_limits_nct_minimize_scalar(\n", + " ncp: float,\n", + " dof: int,\n", + " conf_level: Optional[float] = None,\n", + " alpha_lower: Optional[float] = None,\n", + " alpha_upper: Optional[float] = None,\n", + " tol: float = 1e-9,\n", + " max_iter: int = 10000,\n", + "):\n", + " \"\"\"Calculate confidence limits for the noncentrality parameter (NCP) of the noncentral t-distribution using scipy.optimize.minimize_scalar.\n", + "\n", + " This function uses the scipy.optimize.minimize_scalar method to estimate\n", + " the lower and upper confidence limits for the noncentrality\n", + " parameter (NCP) given the degrees of freedom, confidence level, and other parameters.\n", + "\n", + " Args:\n", + " ncp (float, optional): The observed noncentrality parameter. Can be passed as 't_value'.\n", + " df (int): Degrees of freedom. Must be positive.\n", + " conf_level (float, optional): The confidence level for the interval.\n", + " alpha_lower (float, optional): The significance level for the lower tail.\n", + " alpha_upper (float, optional): The significance level for the upper tail.\n", + " tol (float, optional): Tolerance for the optimization algorithms. Default is 1e-9.\n", + " max_iter (int, optional): Maximum number of iterations to perform. Default is 10000.\n", + "\n", + " Returns:\n", + " dict: A dictionary with the following keys:\n", + " - lower_limit (float): Lower confidence limit for the NCP.\n", + " - prob_less_lower (float): Probability that the NCP is less than the lower limit.\n", + " - upper_limit (float): Upper confidence limit for the NCP.\n", + " - prob_greater_upper (float): Probability that the NCP is greater than the upper limit.\n", + "\n", + " Raises:\n", + " ValueError:\n", "\n", + " Example:\n", + " >>> conf_limits_nct_minimize_scalar(ncp=2.83, dof=126, conf_level=0.95)\n", + " {'lower_limit': 0.8337502600175457, 'prob_less_lower': 0.02499999995262825, 'upper_limit': 4.815359140504376, 'prob_greater_upper': 0.024999999971943743}\n", + " \"\"\"\n", + " # If conf_level is provided without alpha values, calculate them\n", + " if conf_level is not None and alpha_lower is None and alpha_upper is None:\n", + " alpha_lower = (1 - conf_level) / 2\n", + " alpha_upper = (1 - conf_level) / 2\n", + "\n", + " # Internal function to compute lower confidence limit\n", + " def ci_nct_lower(val_of_interest):\n", + " return (nct.ppf(1 - alpha_lower, dof, val_of_interest, loc=0) - ncp) ** 2\n", + "\n", + " # Internal function to compute upper confidence limit\n", + " def ci_nct_upper(val_of_interest):\n", + " return (nct.ppf(alpha_upper, dof, val_of_interest, loc=0) - ncp) ** 2\n", "\n", + " min_ncp = min(-150, -5 * ncp)\n", + " max_ncp = max(150, 5 * ncp)\n", + "\n", + " lower_limit = minimize_scalar(\n", + " ci_nct_lower,\n", + " bounds=(min_ncp, max_ncp),\n", + " method=\"bounded\",\n", + " options={\"xatol\": tol},\n", + " )\n", + " upper_limit = minimize_scalar(\n", + " ci_nct_upper,\n", + " bounds=(min_ncp, max_ncp),\n", + " method=\"bounded\",\n", + " options={\"xatol\": tol, \"disp\": 0, \"maxiter\": max_iter},\n", + " )\n", + "\n", + " return {\n", + " \"lower_limit\": lower_limit.x if alpha_lower != 0 else -np.inf,\n", + " \"prob_less_lower\": (\n", + " 1 - nct.cdf(ncp, dof, lower_limit.x, loc=0) if alpha_lower != 0 else 0\n", + " ),\n", + " \"upper_limit\": upper_limit.x if alpha_upper != 0 else np.inf,\n", + " \"prob_greater_upper\": (\n", + " nct.cdf(ncp, dof, upper_limit.x, loc=0) if alpha_upper != 0 else 0\n", + " ),\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'lower_limit': 0.8337502600175457,\n", + " 'prob_less_lower': 0.02499999995262825,\n", + " 'upper_limit': 4.815359140504376,\n", + " 'prob_greater_upper': 0.024999999971943743}" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "conf_limits_nct_minimize_scalar(ncp=2.83, dof=126, conf_level=0.95)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'lower_limit': 0.46169197879015583,\n", + " 'prob_less_lower': 0.01000000006700108,\n", + " 'upper_limit': 4.602743339958727,\n", + " 'prob_greater_upper': 0.03999999998737774}" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "conf_limits_nct_minimize_scalar(\n", + " ncp=2.83, dof=126, conf_level=None, alpha_lower=0.01, alpha_upper=0.04\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "c:\\projects\\pycvcqv\\.venv\\Lib\\site-packages\\scipy\\optimize\\_optimize.py:2289: RuntimeWarning: invalid value encountered in scalar subtract\n", + " r = (xf - nfc) * (fx - ffulc)\n", + "c:\\projects\\pycvcqv\\.venv\\Lib\\site-packages\\scipy\\optimize\\_optimize.py:2290: RuntimeWarning: invalid value encountered in scalar subtract\n", + " q = (xf - fulc) * (fx - fnfc)\n" + ] + }, + { + "data": { + "text/plain": [ + "{'lower_limit': -inf,\n", + " 'prob_less_lower': 0,\n", + " 'upper_limit': 4.495224841782837,\n", + " 'prob_greater_upper': 0.04999999999219518}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "conf_limits_nct_minimize_scalar(\n", + " ncp=2.83, dof=126, conf_level=None, alpha_lower=0, alpha_upper=0.05\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 133, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "False" + ] + }, + "execution_count": 133, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "0.9750000000280563 > 0.9750000000473718" + ] + }, + { + "cell_type": "code", + "execution_count": 103, + "metadata": {}, + "outputs": [], + "source": [ "def conf_limits_nct(\n", - " ncp: float = None,\n", - " df: int = None,\n", - " conf_level: float = 0.95,\n", - " alpha_lower: float = None,\n", - " alpha_upper: float = None,\n", - " t_value: float = None,\n", + " ncp: float,\n", + " df: int,\n", + " conf_level: Optional[float] = None,\n", + " alpha_lower: Optional[float] = None,\n", + " alpha_upper: Optional[float] = None,\n", " tol: float = 1e-9,\n", - " sup_int_warns: bool = True,\n", ") -> dict:\n", " \"\"\"Calculate confidence limits for the noncentrality parameter (NCP) of the noncentral t-distribution.\n", "\n", @@ -71,9 +390,7 @@ " conf_level (float, optional): The confidence level for the interval. Default is 0.95.\n", " alpha_lower (float, optional): The significance level for the lower tail.\n", " alpha_upper (float, optional): The significance level for the upper tail.\n", - " t_value (float, optional): Alias for 'ncp'. Either 'ncp' or 't_value' must be specified.\n", " tol (float, optional): Tolerance for the optimization algorithms. Default is 1e-9.\n", - " sup_int_warns (bool, optional): Whether to suppress warnings during optimization. Default is True.\n", "\n", " Returns:\n", " dict: A dictionary with the following keys:\n", @@ -90,22 +407,6 @@ " {'lower_limit': 1.97, 'prob_less_lower': 0.025, 'upper_limit': 3.67, 'prob_greater_upper': 0.975}\n", " \"\"\"\n", "\n", - " # If ncp is not provided, use t_value as its alias\n", - " if ncp is None:\n", - " if t_value is None:\n", - " raise ValueError(\"You need to specify either 'ncp' or its alias 't_value'\")\n", - " ncp = t_value\n", - "\n", - " # Check if degrees of freedom are positive\n", - " if df <= 0:\n", - " raise ValueError(\"The degrees of freedom must be some positive value.\")\n", - "\n", - " # Warning if the noncentrality parameter is beyond the accuracy limit\n", - " if abs(ncp) > 37.62:\n", - " print(\n", - " \"Warning: The observed noncentrality parameter exceeds 37.62, which may affect accuracy.\"\n", - " )\n", - "\n", " # If conf_level is provided without alpha values, calculate them\n", " if conf_level is not None and alpha_lower is None and alpha_upper is None:\n", " alpha_lower = (1 - conf_level) / 2\n", @@ -125,10 +426,16 @@ " max_ncp = max(150, 5 * ncp)\n", "\n", " lower_limit = minimize_scalar(\n", - " ci_nct_lower, bounds=(min_ncp, max_ncp), method=\"bounded\", tol=tol\n", + " ci_nct_lower,\n", + " bounds=(min_ncp, max_ncp),\n", + " method=\"bounded\",\n", + " options={\"xatol\": tol},\n", " )\n", " upper_limit = minimize_scalar(\n", - " ci_nct_upper, bounds=(min_ncp, max_ncp), method=\"bounded\", tol=tol\n", + " ci_nct_upper,\n", + " bounds=(min_ncp, max_ncp),\n", + " method=\"bounded\",\n", + " options={\"xatol\": tol},\n", " )\n", "\n", " return {\n", @@ -201,6 +508,29 @@ " }" ] }, + { + "cell_type": "code", + "execution_count": 104, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'lower_limit': 4.815359140504376,\n", + " 'prob_less_lower': 0.024999999971943743,\n", + " 'upper_limit': 0.8337502600175457,\n", + " 'prob_greater_upper': 0.02499999995262825}" + ] + }, + "execution_count": 104, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "conf_limits_nct(ncp=2.83, df=126, conf_level=0.95)" + ] + }, { "cell_type": "code", "execution_count": 3, diff --git a/poetry.lock b/poetry.lock index 5a79c9f..0f70cdd 100644 --- a/poetry.lock +++ b/poetry.lock @@ -3084,6 +3084,48 @@ setuptools = ">=19.3" github = ["jinja2 (>=3.1.0)", "pygithub (>=1.43.3)"] gitlab = ["python-gitlab (>=1.3.0)"] +[[package]] +name = "scipy" +version = "1.13.1" +description = "Fundamental algorithms for scientific computing in Python" +optional = false +python-versions = ">=3.9" +files = [ + {file = "scipy-1.13.1-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:20335853b85e9a49ff7572ab453794298bcf0354d8068c5f6775a0eabf350aca"}, + {file = "scipy-1.13.1-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:d605e9c23906d1994f55ace80e0125c587f96c020037ea6aa98d01b4bd2e222f"}, + {file = "scipy-1.13.1-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:cfa31f1def5c819b19ecc3a8b52d28ffdcc7ed52bb20c9a7589669dd3c250989"}, + {file = "scipy-1.13.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:f26264b282b9da0952a024ae34710c2aff7d27480ee91a2e82b7b7073c24722f"}, + {file = "scipy-1.13.1-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:eccfa1906eacc02de42d70ef4aecea45415f5be17e72b61bafcfd329bdc52e94"}, + {file = "scipy-1.13.1-cp310-cp310-win_amd64.whl", hash = "sha256:2831f0dc9c5ea9edd6e51e6e769b655f08ec6db6e2e10f86ef39bd32eb11da54"}, + {file = "scipy-1.13.1-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:27e52b09c0d3a1d5b63e1105f24177e544a222b43611aaf5bc44d4a0979e32f9"}, + {file = "scipy-1.13.1-cp311-cp311-macosx_12_0_arm64.whl", hash = "sha256:54f430b00f0133e2224c3ba42b805bfd0086fe488835effa33fa291561932326"}, + {file = "scipy-1.13.1-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:e89369d27f9e7b0884ae559a3a956e77c02114cc60a6058b4e5011572eea9299"}, + {file = "scipy-1.13.1-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:a78b4b3345f1b6f68a763c6e25c0c9a23a9fd0f39f5f3d200efe8feda560a5fa"}, + {file = "scipy-1.13.1-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:45484bee6d65633752c490404513b9ef02475b4284c4cfab0ef946def50b3f59"}, + {file = "scipy-1.13.1-cp311-cp311-win_amd64.whl", hash = "sha256:5713f62f781eebd8d597eb3f88b8bf9274e79eeabf63afb4a737abc6c84ad37b"}, + {file = "scipy-1.13.1-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:5d72782f39716b2b3509cd7c33cdc08c96f2f4d2b06d51e52fb45a19ca0c86a1"}, + {file = "scipy-1.13.1-cp312-cp312-macosx_12_0_arm64.whl", hash = "sha256:017367484ce5498445aade74b1d5ab377acdc65e27095155e448c88497755a5d"}, + {file = "scipy-1.13.1-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:949ae67db5fa78a86e8fa644b9a6b07252f449dcf74247108c50e1d20d2b4627"}, + {file = "scipy-1.13.1-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:de3ade0e53bc1f21358aa74ff4830235d716211d7d077e340c7349bc3542e884"}, + {file = "scipy-1.13.1-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:2ac65fb503dad64218c228e2dc2d0a0193f7904747db43014645ae139c8fad16"}, + {file = "scipy-1.13.1-cp312-cp312-win_amd64.whl", hash = "sha256:cdd7dacfb95fea358916410ec61bbc20440f7860333aee6d882bb8046264e949"}, + {file = "scipy-1.13.1-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:436bbb42a94a8aeef855d755ce5a465479c721e9d684de76bf61a62e7c2b81d5"}, + {file = "scipy-1.13.1-cp39-cp39-macosx_12_0_arm64.whl", hash = "sha256:8335549ebbca860c52bf3d02f80784e91a004b71b059e3eea9678ba994796a24"}, + {file = "scipy-1.13.1-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:d533654b7d221a6a97304ab63c41c96473ff04459e404b83275b60aa8f4b7004"}, + {file = "scipy-1.13.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:637e98dcf185ba7f8e663e122ebf908c4702420477ae52a04f9908707456ba4d"}, + {file = "scipy-1.13.1-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:a014c2b3697bde71724244f63de2476925596c24285c7a637364761f8710891c"}, + {file = "scipy-1.13.1-cp39-cp39-win_amd64.whl", hash = "sha256:392e4ec766654852c25ebad4f64e4e584cf19820b980bc04960bca0b0cd6eaa2"}, + {file = "scipy-1.13.1.tar.gz", hash = "sha256:095a87a0312b08dfd6a6155cbbd310a8c51800fc931b8c0b84003014b874ed3c"}, +] + +[package.dependencies] +numpy = ">=1.22.4,<2.3" + +[package.extras] +dev = ["cython-lint (>=0.12.2)", "doit (>=0.36.0)", "mypy", "pycodestyle", "pydevtool", "rich-click", "ruff", "types-psutil", "typing_extensions"] +doc = ["jupyterlite-pyodide-kernel", "jupyterlite-sphinx (>=0.12.0)", "jupytext", "matplotlib (>=3.5)", "myst-nb", "numpydoc", "pooch", "pydata-sphinx-theme (>=0.15.2)", "sphinx (>=5.0.0)", "sphinx-design (>=0.4.0)"] +test = ["array-api-strict", "asv", "gmpy2", "hypothesis (>=6.30)", "mpmath", "pooch", "pytest", "pytest-cov", "pytest-timeout", "pytest-xdist", "scikit-umfpack", "threadpoolctl"] + [[package]] name = "setuptools" version = "72.2.0" @@ -3880,4 +3922,4 @@ test = ["big-O", "importlib-resources", "jaraco.functools", "jaraco.itertools", [metadata] lock-version = "2.0" python-versions = ">=3.9,<4.0" -content-hash = "1141ef4452318b36993f4d2d5ed716c00c40053f41754b2b5db8999880d488b1" +content-hash = "05875d1d1b3e362adba7a3b028fcbcc2405105c10eaebd56657eddda4757b727" diff --git a/pycvcqv/noncentralt.py b/pycvcqv/noncentralt.py new file mode 100644 index 0000000..803b432 --- /dev/null +++ b/pycvcqv/noncentralt.py @@ -0,0 +1,112 @@ +"""Noncentral t-distribution module.""" + +# --------------------------- Import libraries and functions -------------------------- +from typing import Dict, Optional, Union + +import numpy as np +from scipy.optimize import minimize_scalar +from scipy.stats import nct + +from pycvcqv.checkers import is_dof_positive_natural_number, is_ncp_huge +from pycvcqv.sanitizers import validate_ncp_confidence_level_arguments + + +# -------------------------------- function definition -------------------------------- +@is_dof_positive_natural_number +@is_ncp_huge +@validate_ncp_confidence_level_arguments +def conf_limits_nct_minimize_scalar( + ncp: float, + dof: int, + conf_level: Optional[float] = None, + alpha_lower: Optional[float] = None, + alpha_upper: Optional[float] = None, + tol: Optional[float] = 1e-9, + max_iter: Optional[int] = 10000, +) -> Dict[str, Union[float, int]]: + """ + Calculate confidence limits for the noncentrality parameter (NCP) of the + noncentral t-distribution using scipy.optimize.minimize_scalar. + + This function uses the scipy.optimize.minimize_scalar method to estimate + the lower and upper confidence limits for the noncentrality + parameter (NCP) given the degrees of freedom, confidence level, and other parameters. + + Args: + ncp (float): The observed noncentrality parameter. Can be passed as 't_value'. + dof (int): Degrees of freedom. Must be positive. + conf_level (float, optional): The confidence level for the interval. + alpha_lower (float, optional): The significance level for the lower tail. + alpha_upper (float, optional): The significance level for the upper tail. + tol (float, optional): Tolerance for the optimization algorithms. Default is 1e-9. + max_iter (int, optional): Maximum number of iterations to perform. Default is 10000. + + Returns: + dict: A dictionary with the following keys: + - lower_limit (float): Lower confidence limit for the NCP. + - prob_less_lower (float): Probability that the NCP is less than the lower limit. + - upper_limit (float): Upper confidence limit for the NCP. + - prob_greater_upper (float): Probability that the NCP is greater than the upper limit. + + Example: + .. code:: python + >>> conf_limits_nct_minimize_scalar(ncp=2.83, dof=126, conf_level=0.95) + ... { + ... 'lower_limit': 0.8337502600175457, + ... 'prob_less_lower': 0.02499999995262825, + ... 'upper_limit': 4.815359140504376, + ... 'prob_greater_upper': 0.024999999971943743 + ... } + """ + # --- If all three are None, use default conf_level and compute alpha values -- + if conf_level is None and alpha_lower is None and alpha_upper is None: + conf_level = 0.95 + alpha_lower = (1 - conf_level) / 2 + alpha_upper = (1 - conf_level) / 2 + # ---- Calculate the alpha_lower and alpha_upper based on given conf_level ---- + elif conf_level is not None and alpha_lower is None and alpha_upper is None: + alpha_lower = (1 - conf_level) / 2 + alpha_upper = (1 - conf_level) / 2 + + def _ci_nct_lower(val_of_interest: float) -> float: + """Internal function to compute lower confidence limit.""" + assert alpha_lower is not None # Ensuring alpha_lower is not None + result: float = nct.ppf( + 1 - alpha_lower, dof, val_of_interest, loc=0 + ) # Explicit type declaration + return (result - ncp) ** 2 + + def _ci_nct_upper(val_of_interest: float) -> float: + """Internal function to compute upper confidence limit.""" + assert alpha_upper is not None # Ensuring alpha_upper is not None + result: float = nct.ppf( + alpha_upper, dof, val_of_interest, loc=0 + ) # Explicit type declaration + return (result - ncp) ** 2 + + min_ncp = min(-150, -5 * ncp) + max_ncp = max(150, 5 * ncp) + + lower_limit = minimize_scalar( + _ci_nct_lower, + bounds=(min_ncp, max_ncp), + method="bounded", + options={"xatol": tol}, + ) + upper_limit = minimize_scalar( + _ci_nct_upper, + bounds=(min_ncp, max_ncp), + method="bounded", + options={"xatol": tol, "disp": 0, "maxiter": max_iter}, + ) + + return { + "lower_limit": lower_limit.x if alpha_lower != 0 else -np.inf, + "prob_less_lower": ( + 1 - nct.cdf(ncp, dof, lower_limit.x, loc=0) if alpha_lower != 0 else 0 + ), + "upper_limit": upper_limit.x if alpha_upper != 0 else np.inf, + "prob_greater_upper": ( + nct.cdf(ncp, dof, upper_limit.x, loc=0) if alpha_upper != 0 else 0 + ), + } diff --git a/pycvcqv/sanitizers.py b/pycvcqv/sanitizers.py index 9d2e239..11bc845 100644 --- a/pycvcqv/sanitizers.py +++ b/pycvcqv/sanitizers.py @@ -15,17 +15,6 @@ def wrapper(*args, **kwargs): alpha_lower = kwargs.get("alpha_lower", args[3] if len(args) > 3 else None) alpha_upper = kwargs.get("alpha_upper", args[4] if len(args) > 4 else None) - # --- If all three are None, use default conf_level and compute alpha values -- - if conf_level is None and alpha_lower is None and alpha_upper is None: - conf_level = 0.95 - alpha_lower = (1 - conf_level) / 2 - alpha_upper = (1 - conf_level) / 2 - - # ---- Calculate the alpha_lower and alpha_upper based on given conf_level ---- - elif conf_level is not None and alpha_lower is None and alpha_upper is None: - alpha_lower = (1 - conf_level) / 2 - alpha_upper = (1 - conf_level) / 2 - # -- Validate if one of alpha_lower or alpha_upper is given, both must be given if (alpha_lower is not None) != (alpha_upper is not None): raise ValueError( @@ -33,11 +22,21 @@ def wrapper(*args, **kwargs): ) # --- Validate that alpha_lower and alpha_upper are within the [0, 1] range --- - if alpha_lower is not None and not 0 <= alpha_lower <= 1: - raise ValueError("conf_level or alpha values must be in the range [0, 1].") - if alpha_upper is not None and not 0 <= alpha_upper <= 1: - raise ValueError("conf_level or alpha values must be in the range [0, 1].") - + if alpha_lower is not None: + if not 0 <= alpha_lower <= 1: + raise ValueError( + "conf_level and alpha values must be in the range [0, 1]." + ) + if alpha_upper is not None: + if not 0 <= alpha_upper <= 1: + raise ValueError( + "conf_level and alpha values must be in the range [0, 1]." + ) + if conf_level is not None: + if not 0 <= conf_level <= 1: + raise ValueError( + "conf_level and alpha values must be in the range [0, 1]." + ) # --------------- Update the kwargs with validated alpha values --------------- kwargs.update( { diff --git a/pyproject.toml b/pyproject.toml index a029be8..3e19ff8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,7 @@ classifiers = [ numpy = "^1.26.4" pandas = "*" python = ">=3.9,<4.0" +scipy = "*" [tool.poetry.dev-dependencies] bandit = "^1.7.4" diff --git a/requirements.txt b/requirements.txt index 1f29e9c6f3f978ef6fd7f2d92fb4d367507e7be5..2e7eaaf3fada217ed959a76ab7ef99da19adf02e 100644 GIT binary patch delta 32 kcmX@hc9ngD1CwwvLncE3LnVVP5E?S*0kQGs(~M1w0EvYMsQ>@~ delta 12 Tcmcc0ewJ;61JmXcOuUQ$AejU$ diff --git a/tests/test_checkers.py b/tests/test_checkers.py index 75afe0b..c5d47ba 100644 --- a/tests/test_checkers.py +++ b/tests/test_checkers.py @@ -111,13 +111,6 @@ def test_default_conf_level(): assert result == (2.0 * 100) * (expected_alpha * 2) -def test_given_conf_level(): - """Test case where conf_level is provided, but alphas are not.""" - result = mocked_conf_limits_nct(2.0, 100, conf_level=0.90) - expected_alpha = (1 - 0.90) / 2 - assert result == (2.0 * 100) * (expected_alpha * 2) - - def test_given_alpha_values(): """Test case where alpha_lower and alpha_upper are provided.""" result = mocked_conf_limits_nct(2.0, 100, alpha_lower=0.01, alpha_upper=0.02) @@ -140,7 +133,7 @@ def test_invalid_alpha_lower(): mocked_conf_limits_nct(2.0, 100, alpha_lower=-0.1, alpha_upper=0.02) assert ( execinfo.value.args[0] - == "conf_level or alpha values must be in the range [0, 1]." + == "conf_level and alpha values must be in the range [0, 1]." ) @@ -150,15 +143,17 @@ def test_invalid_alpha_upper(): mocked_conf_limits_nct(2.0, 100, alpha_lower=0.01, alpha_upper=1.2) assert ( execinfo.value.args[0] - == "conf_level or alpha values must be in the range [0, 1]." + == "conf_level and alpha values must be in the range [0, 1]." ) def test_invalid_conf_level(): """Test case where conf_level is out of the [0, 1] range.""" with pytest.raises(ValueError) as execinfo: - mocked_conf_limits_nct(2.0, 100, conf_level=1.5) + mocked_conf_limits_nct( + ncp=2.0, dof=100, conf_level=1.5, alpha_lower=0.01, alpha_upper=0.95 + ) assert ( execinfo.value.args[0] - == "conf_level or alpha values must be in the range [0, 1]." + == "conf_level and alpha values must be in the range [0, 1]." ) diff --git a/tests/test_noncentralt.py b/tests/test_noncentralt.py new file mode 100644 index 0000000..0f25ac2 --- /dev/null +++ b/tests/test_noncentralt.py @@ -0,0 +1,53 @@ +"""Tests for noncentral t-distribution module.""" + +# --------------------------- Import libraries and functions -------------------------- + +from math import isinf + +import pytest + +from pycvcqv.noncentralt import conf_limits_nct_minimize_scalar + + +def test_conf_limits_nct_minimize_scalar_all_none(): + """Tests the conf_limits_nct_minimize_scalar when conf_level and alphas are None.""" + result = conf_limits_nct_minimize_scalar(ncp=2.83, dof=126) + + assert result["lower_limit"] == pytest.approx(0.8337502600175457, rel=1e-9) + assert result["prob_less_lower"] == pytest.approx(0.02499999995262825, rel=1e-9) + assert result["upper_limit"] == pytest.approx(4.815359140504376, rel=1e-9) + assert result["prob_greater_upper"] == pytest.approx(0.024999999971943743, rel=1e-9) + + +def test_conf_limits_nct_minimize_scalar_conf_level_95(): + """Tests the conf_limits_nct_minimize_scalar when conf_level is set to 0.95.""" + result = conf_limits_nct_minimize_scalar(ncp=2.83, dof=126, conf_level=0.95) + + assert result["lower_limit"] == pytest.approx(0.8337502600175457, rel=1e-9) + assert result["prob_less_lower"] == pytest.approx(0.02499999995262825, rel=1e-9) + assert result["upper_limit"] == pytest.approx(4.815359140504376, rel=1e-9) + assert result["prob_greater_upper"] == pytest.approx(0.024999999971943743, rel=1e-9) + + +def test_conf_limits_nct_minimize_scalar_alpha_lower_0_01_alpha_upper_0_04(): + """Tests the conf_limits_nct_minimize_scalar when alpha levels are non-zero.""" + result = conf_limits_nct_minimize_scalar( + ncp=2.83, dof=126, conf_level=None, alpha_lower=0.01, alpha_upper=0.04 + ) + + assert result["lower_limit"] == pytest.approx(0.46169197879015583, rel=1e-9) + assert result["prob_less_lower"] == pytest.approx(0.01000000006700108, rel=1e-9) + assert result["upper_limit"] == pytest.approx(4.602743339958727, rel=1e-9) + assert result["prob_greater_upper"] == pytest.approx(0.03999999998737774, rel=1e-9) + + +def test_conf_limits_nct_minimize_scalar_alpha_lower_0_alpha_upper_0_05(): + """Tests the conf_limits_nct_minimize_scalar when alpha_lower is zero.""" + result = conf_limits_nct_minimize_scalar( + ncp=2.83, dof=126, conf_level=None, alpha_lower=0, alpha_upper=0.05 + ) + + assert isinf(result["lower_limit"]) and result["lower_limit"] < 0 + assert result["prob_less_lower"] == 0 + assert result["upper_limit"] == pytest.approx(4.495224841782837, rel=1e-9) + assert result["prob_greater_upper"] == pytest.approx(0.04999999999219518, rel=1e-9)