From cb2e3a598138aca5a55ed67c4d77cc20c72a80b4 Mon Sep 17 00:00:00 2001 From: saudzahirr Date: Sat, 15 Aug 2026 14:07:52 +0500 Subject: [PATCH 1/2] Fix the second-order moment equations to match Cox (1979) rawmoment() diverged from the equations in the Appendix of N. D. Cox, "Tolerance Analysis by Computer", Journal of Quality Technology 11(2), 1979. The mean and variance (A-6, A-7) were correct, but the third raw moment was wrong for n >= 3 inputs and the fourth for n >= 2. Nine defects, each checked against the printed equations on pp. 86-87 and against two independent oracles (exact polynomial expansion of E[y^k], and tensor-product Gauss-Hermite/Gauss-Laguerre quadrature): 1. (A-8) triple sum used vm[m,1] instead of vm[m,2]. The first standardized moment is zero, so this silently dropped the term. 2. (A-9) pair term used mu_j6 instead of mu_j5 in 12 b_ij b_ii b_jj [b_ii mu_i5 mu_j3 + b_jj mu_i3 mu_j5]. 3-5. (A-9) three terms indexed mu_i2 where the paper has mu_j2. 6. (A-9) triple sum used 6 b_ii^2 b_jk instead of 6 b_ii^2 b_jk^2. 7. (A-9) triple sum used 24 b_ii b_kk b_ik b_jk instead of 24 b_ii b_jj b_ik b_jk. 8. (A-9) emitted two mu_i2 mu_j3 mu_k3 blocks, double counting three terms; the paper has one. 9. (A-9) quadruple sum used b_jj b_ik b_im b_im instead of b_jj b_ik b_im b_km. With these corrected, Cox's Table 2 reproduces: mean 1176.45, variance 99,699.68, and skewness 0.7100769 against the paper's 0.7100789 (it was 0.7080131). Table 2's printed kurtosis of 6.183188 is not reproducible from equation (A-9) itself, which gives 6.161161; the code follows the equation, and test_assembly.py records why. Also fixed, found while verifying the surrounding layers: - umath.cot() computed tan() where it needed csc(), so both the first and second derivatives were wrong. - UncertainVariable(rv=...) ignored loc/scale and skipped standardization for distributions with no shape parameter, so U(a, b) returned unstandardized moments (mu_4 = 0.0125 rather than 1.8). Both branches now integrate the standardized central moments directly, which also avoids the cancellation the raw-moment route suffers when the mean is large relative to the spread. docs/theory.md carried five equation errors of its own, independent of the code; these are corrected, and the blocks previously elided as "similar terms" are written out in full. Stale published results in README.md, soerp_examples.py and the uv docstring are updated. Co-authored-by: codex <267193182+codex@users.noreply.github.com> Co-authored-by: Claude Opus 5 Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- README.md | 8 ++--- docs/theory.md | 60 ++++++++++++++++++++++++++++-------- soerp/method_of_moments.py | 27 +++++----------- soerp/umath.py | 4 +-- soerp/uncertain_variable.py | 61 +++++++++++++++++-------------------- soerp_examples.py | 12 ++++---- tests/test_assembly.py | 9 ++++-- 7 files changed, 102 insertions(+), 79 deletions(-) diff --git a/README.md b/README.md index 29e6a95..87f45f5 100644 --- a/README.md +++ b/README.md @@ -64,7 +64,7 @@ We can now see the results of the calculations in two ways: ```python >>> Z # "print" is optional at the command-line -uv(1176.45, 99699.6822917, 0.708013052944, 6.16324345127) +uv(1176.45, 99699.6822917, 0.710076903845, 6.16116115559) ``` 2. The ``describe`` class method that explains briefly what the values are: @@ -74,8 +74,8 @@ uv(1176.45, 99699.6822917, 0.708013052944, 6.16324345127) SOERP Uncertain Value: > Mean................... 1176.45 > Variance............... 99699.6822917 - > Skewness Coefficient... 0.708013052944 - > Kurtosis Coefficient... 6.16324345127 + > Skewness Coefficient... 0.710076903845 + > Kurtosis Coefficient... 6.16116115559 ``` ### Distribution Moments @@ -86,7 +86,7 @@ The eight moments of any input variable (and four of any output variable) can be >>> x1.moments() [24.0, 1.0, 0.0, 3.0000000000000053, 0.0, 15.000000000000004, 0.0, 105.0] >>> Z.moments() -[1176.45, 99699.6822917, 0.708013052944, 6.16324345127] +[1176.45, 99699.6822917, 0.710076903845, 6.16116115559] ``` ### Correlations diff --git a/docs/theory.md b/docs/theory.md index 02ca9d3..4c8cbf9 100644 --- a/docs/theory.md +++ b/docs/theory.md @@ -170,12 +170,12 @@ $$ \begin{align} \nu_{y3} &= \sum_{i=1}^n \left[b_i^3\,\mu_{i3} + b_{ii}^3\,\mu_{i6} + 3b_i^2 b_{ii}\,\mu_{i4} + 3b_i b_{ii}^2\,\mu_{i5}\right] \\ -&\quad + \sum_{\substack{i,j \\ i < j}} \left[b_{ij}^2\,\mu_{i3}\,\mu_{j3} +&\quad + \sum_{\substack{i,j \\ i < j}} \left[b_{ij}^3\,\mu_{i3}\,\mu_{j3} + 6b_i b_j b_{ij}\,\mu_{i2}\,\mu_{j2} + 6b_{ii}b_{jj}b_{ij}\,\mu_{i3}\,\mu_{j3}\right] \\ &\quad + \sum_{i=1}^n \sum_{\substack{j=1 \\ j \ne i}}^n \left[3b_{ii}^2\,\mu_{i4}\,b_{jj}\,\mu_{j2} - + 6b_i b_{ii} b_{ij}\,\mu_{i2}\,\mu_{j3} + + 6b_i b_{jj} b_{ij}\,\mu_{i2}\,\mu_{j3} + 3b_{ii}b_j^2\,\mu_{i2}\,\mu_{j2} + 6b_i b_{ii} b_{jj}\,\mu_{i3}\,\mu_{j2} + 3b_i b_{ij}^2\,\mu_{i3}\,\mu_{j2} @@ -230,9 +230,16 @@ $$ $$ \sum_{i=1}^n \sum_{\substack{j=1 \\ j \ne i}}^n -\left\{4b_{ii}^3 b_{jj}\,\mu_{i6}\,\mu_{j2} -+ \cdots -+ 6b_i^2 b_{ij}^2\,\mu_{i4}\,\mu_{j2}\right\} +\left\{ + 4b_{ii}^3 b_{jj}\,\mu_{i6}\,\mu_{j2} + + 4b_{ii} b_j^3\,\mu_{i2}\,\mu_{j3} + + 12b_i b_{ii} b_j^2\,\mu_{i3}\,\mu_{j2} + + 12b_i b_{ii}^2 b_{jj}\,\mu_{i5}\,\mu_{j2} + + 12b_i b_{ii} b_{jj}^2\,\mu_{i3}\,\mu_{j4} + + 4b_i b_{ij}^3\,\mu_{i4}\,\mu_{j3} + + 4b_{ii} b_{ij}^3\,\mu_{i5}\,\mu_{j3} + + 6b_{ii}^2 b_j^2\,\mu_{i4}\,\mu_{j2} +\right\} $$ **Triplet terms** ($i < j < k$): @@ -241,21 +248,48 @@ $$ \sum_{i=1}^{n-2}\sum_{j=i+1}^{n-1}\sum_{k=j+1}^{n} \Bigl\{ \left[12b_{ii}^2 b_{jj}b_{kk} + 6b_{ij}^2 b_{ik}^2 - + 12b_{ii}(b_{kk}b_{ij}^2 + b_{jj}b_{ik}^2) + 6b_{ii}^2 b_{jk}^2\right]\mu_{i4}\,\mu_{j2}\,\mu_{k2} - + \left[\text{similar terms for } (j,i,k)\text{ and }(k,i,j)\right] + + 12b_{ii}(b_{kk}b_{ij}^2 + b_{jj}b_{ik}^2) + 6b_{ii}^2 b_{jk}^2\right]\mu_{i4}\,\mu_{j2}\,\mu_{k2} \\ + + \left[12b_{ii} b_{jj}^2 b_{kk} + 6b_{ij}^2 b_{jk}^2 + + 12b_{jj}(b_{kk}b_{ij}^2 + b_{ii}b_{jk}^2) + 6b_{jj}^2 b_{ik}^2\right]\mu_{i2}\,\mu_{j4}\,\mu_{k2} \\ + + \left[12b_{ii} b_{jj} b_{kk}^2 + 6b_{ik}^2 b_{jk}^2 + + 12b_{kk}(b_{ii}b_{jk}^2 + b_{jj}b_{ik}^2) + 6b_{kk}^2 b_{ij}^2\right]\mu_{i2}\,\mu_{j2}\,\mu_{k4} \\ + \left[12b_{ij}^2 b_{ik}b_{jk} + 24b_{ii}b_{jj}b_{kk}b_{ij} - + 4b_{kk}b_{ij}^3 + 24b_{ii}b_{jj}b_{ik}b_{jk}\right]\mu_{i3}\,\mu_{j3}\,\mu_{k2} - + \left[\text{similar terms for remaining } \mu_{i3}\mu_{j2}\mu_{k3}\text{ and }\mu_{i2}\mu_{j3}\mu_{k3}\right] + + 4b_{kk}b_{ij}^3 + 24b_{ii}b_{jj}b_{ik}b_{jk}\right]\mu_{i3}\,\mu_{j3}\,\mu_{k2} \\ + + \left[12b_{ij} b_{ik}^2 b_{jk} + 24b_{ii}b_{jj}b_{kk}b_{ik} + + 4b_{jj}b_{ik}^3 + 24b_{ii}b_{kk}b_{ij}b_{jk}\right]\mu_{i3}\,\mu_{j2}\,\mu_{k3} \\ + + \left[12b_{ij} b_{ik} b_{jk}^2 + 24b_{ii}b_{jj}b_{kk}b_{jk} + + 4b_{ii}b_{jk}^3 + 24b_{jj}b_{kk}b_{ij}b_{ik}\right]\mu_{i2}\,\mu_{j3}\,\mu_{k3} \\ + 24\left[b_{ii}b_{jj}b_{kk} + b_{ij}b_{ik}b_{jk}\right] \cdot\left[b_i\,\mu_{i3}\,\mu_{j2}\,\mu_{k2} + b_j\,\mu_{i2}\,\mu_{j3}\,\mu_{k2} - + b_k\,\mu_{i2}\,\mu_{j2}\,\mu_{k3}\right] + + b_k\,\mu_{i2}\,\mu_{j2}\,\mu_{k3}\right] \\ + + 12\left[b_i b_{jk}^2\,\mu_{i2}(b_{ij}\,\mu_{j3}\,\mu_{k2} + b_{ik}\,\mu_{j2}\,\mu_{k3}) + + b_j b_{ik}^2\,\mu_{j2}(b_{ij}\,\mu_{i3}\,\mu_{k2} + b_{jk}\,\mu_{i2}\,\mu_{k3}) + + b_k b_{ij}^2\,\mu_{k2}(b_{ik}\,\mu_{i3}\,\mu_{j2} + b_{jk}\,\mu_{i2}\,\mu_{j3})\right] \\ + + 12\left[b_{ii} b_{jk}^2\,\mu_{i3}(b_{ij}\,\mu_{j3}\,\mu_{k2} + b_{ik}\,\mu_{j2}\,\mu_{k3}) + + b_{jj} b_{ik}^2\,\mu_{j3}(b_{ij}\,\mu_{i3}\,\mu_{k2} + b_{jk}\,\mu_{i2}\,\mu_{k3}) + + b_{kk} b_{ij}^2\,\mu_{k3}(b_{ik}\,\mu_{i3}\,\mu_{j2} + b_{jk}\,\mu_{i2}\,\mu_{j3})\right] \\ + + 24b_{ij}b_{ik}b_{jk}\left[b_{ii}\,\mu_{i4}\,\mu_{j2}\,\mu_{k2} + + b_{jj}\,\mu_{i2}\,\mu_{j4}\,\mu_{k2} + + b_{kk}\,\mu_{i2}\,\mu_{j2}\,\mu_{k4}\right] \\ + \mu_{i2}\,\mu_{j2}\,\mu_{k2}\Bigl[ - 12\left(b_{ii}b_{jj}b_{kk}^2 + b_{ii}b_{kk}b_{jk}^2 + b_{jj}b_{kk}b_{ik}^2\right) \\ - + 6\left(b_{ii}^2 b_{jk}^2 + b_{jj}^2 b_{ik}^2 + b_{kk}^2 b_{ij}^2\right) + 12\left(b_{ii}b_{jj}b_k^2 + b_{ii}b_{kk}b_j^2 + b_{jj}b_{kk}b_i^2\right) + + 6\left(b_i^2 b_{jk}^2 + b_j^2 b_{ik}^2 + b_k^2 b_{ij}^2\right) \\ + 24\left(b_{ij}b_{ik}b_j b_k + b_{ij}b_{jk}b_i b_k + b_{ik}b_{jk}b_i b_j\right) + 24\left(b_i b_j b_{kk}b_{ij} + b_i b_k b_{jj}b_{ik} + b_j b_k b_{ii}b_{jk}\right) - \Bigr] + \Bigr] \\ + + \mu_{i3}\,\mu_{j2}\,\mu_{k2}\Bigl[ + 24b_j b_{ij} b_{ii} b_{kk} + 24b_k b_{ik} b_{ii} b_{jj} + 12b_i b_{jk}^2 b_{ii} + + 24b_j b_{ik} b_{jk} b_{ii} + 24b_k b_{ij} b_{jk} b_{ii} + + 12b_i b_{ik}^2 b_{jj} + 12b_i b_{ij}^2 b_{kk}\Bigr] \\ + + \mu_{i2}\,\mu_{j3}\,\mu_{k2}\Bigl[ + 24b_i b_{ij} b_{jj} b_{kk} + 24b_k b_{jk} b_{ii} b_{jj} + 12b_j b_{ik}^2 b_{jj} + + 24b_i b_{ik} b_{jk} b_{jj} + 24b_k b_{ij} b_{ik} b_{jj} + + 12b_j b_{jk}^2 b_{ii} + 12b_j b_{ij}^2 b_{kk}\Bigr] \\ + + \mu_{i2}\,\mu_{j2}\,\mu_{k3}\Bigl[ + 24b_i b_{ik} b_{jj} b_{kk} + 24b_j b_{jk} b_{ii} b_{kk} + 12b_k b_{ij}^2 b_{kk} + + 24b_i b_{ij} b_{jk} b_{kk} + 24b_j b_{ij} b_{ik} b_{kk} + + 12b_k b_{jk}^2 b_{ii} + 12b_k b_{ik}^2 b_{jj}\Bigr] \Bigr\} $$ diff --git a/soerp/method_of_moments.py b/soerp/method_of_moments.py index a3f6497..930b4c7 100644 --- a/soerp/method_of_moments.py +++ b/soerp/method_of_moments.py @@ -256,7 +256,7 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] ) * vm[i, 2] * vm[j, 2] - * vm[m, 1] + * vm[m, 2] ) ############################ @@ -291,7 +291,7 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] * qc[j] * ( qc[i] * vm[i, 5] * vm[j, 3] - + qc[j] * vm[i, 3] * vm[j, 6] + + qc[j] * vm[i, 3] * vm[j, 5] ) + 12 * qc[i] @@ -352,7 +352,7 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] * ( qc[i] * ( - lc[i] * vm[i, 5] * vm[i, 2] + lc[i] * vm[i, 5] * vm[j, 2] + lc[j] * vm[i, 4] * vm[j, 3] ) + qc[j] @@ -373,13 +373,13 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] * qc[i] * lc[j] ** 2 * vm[i, 3] - * vm[i, 2] + * vm[j, 2] + 12 * lc[i] * qc[i] ** 2 * qc[j] * vm[i, 5] - * vm[i, 2] + * vm[j, 2] + 12 * lc[i] * qc[i] @@ -405,7 +405,7 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] qc[kk] * cp[i, j] ** 2 + qc[j] * cp[i, kk] ** 2 ) - + 6 * qc[i] ** 2 * cp[j, kk] + + 6 * qc[i] ** 2 * cp[j, kk] ** 2 ) * vm[i, 4] * vm[j, 2] @@ -448,7 +448,7 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] 12 * cp[i, j] ** 2 * cp[i, kk] * cp[j, kk] + 24 * qc[i] * qc[j] * qc[kk] * cp[i, j] + 4 * qc[kk] * cp[i, j] ** 3 - + 24 * qc[i] * qc[kk] * cp[i, kk] * cp[j, kk] + + 24 * qc[i] * qc[j] * cp[i, kk] * cp[j, kk] ) * vm[i, 3] * vm[j, 3] @@ -470,17 +470,6 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] 12 * cp[i, j] * cp[i, kk] * cp[j, kk] ** 2 + 24 * qc[i] * qc[j] * qc[kk] * cp[j, kk] + 4 * qc[i] * cp[j, kk] ** 3 - + 24 * qc[j] * qc[j] * cp[i, j] * cp[i, kk] - ) - * vm[i, 2] - * vm[j, 3] - * vm[kk, 3] - ) - ans += ( - ( - 12 * cp[i, j] * cp[i, kk] * cp[j, kk] - + 24 * qc[i] * qc[j] * qc[kk] * cp[j, kk] - + 4 * qc[i] * cp[j, kk] ** 3 + 24 * qc[j] * qc[kk] * cp[i, j] * cp[i, kk] ) * vm[i, 2] @@ -663,7 +652,7 @@ def rawmoment( # ruff: ignore[too-many-branches, too-many-statements] + qc[j] * cp[i, kk] * cp[i, m] - * cp[i, m] + * cp[kk, m] + qc[kk] * cp[i, j] * cp[i, m] diff --git a/soerp/umath.py b/soerp/umath.py index ffc758e..600d50e 100644 --- a/soerp/umath.py +++ b/soerp/umath.py @@ -138,8 +138,8 @@ def cot(x: UncertainFunction | float) -> UncertainFunction | float: fx = x.x if isinstance(x, UncertainFunction) else x c = math.cos(fx) s = math.sin(fx) - ct = c / s - cs = s / c + ct = c / s # cot(x) + cs = 1.0 / s # csc(x) d1 = -(cs * cs) # $-\csc^2$ d2 = 2.0 * ct * cs * cs # $2 \cot \csc^2$ return _uf_unary(x, ct, d1, d2) diff --git a/soerp/uncertain_variable.py b/soerp/uncertain_variable.py index cae8b0b..87ded88 100644 --- a/soerp/uncertain_variable.py +++ b/soerp/uncertain_variable.py @@ -1,7 +1,6 @@ import numpy as np import scipy.stats as ss -from .method_of_moments import raw2central from .uncertain_function import UncertainFunction @@ -182,7 +181,7 @@ class UncertainVariable(UncertainFunction): >>> x3 = uv([0.5, 0.25, 2, 9, 44, 265, 1854, 14833]) # exp. distributed >>> Z = (x1*x2**2)/(15*(1.5 + x3)) >>> Z - uv(1176.45, 99699.6822919, 0.708013052954, 6.16324345122) + uv(1176.45, 99699.6822919, 0.710076903845, 6.16116115559) The result shows the mean, variance, and standardized skewness and kurtosis of the output variable Z. @@ -273,41 +272,37 @@ def __init__( scale = rv.kwds.get("scale", 1.0) shape = rv.args + if shape and rv.dist.numargs < 1: + raise ValueError( + "The distribution provided doesn't support" + " a 'shape' parameter" + ) + if not shape and rv.dist.numargs != 0: + raise ValueError( + "The distribution provided requires a third" + " 'shape' parameter" + ) + mn = rv.mean() sd = rv.std() - if shape: - if rv.dist.numargs < 1: - raise ValueError( - "The distribution provided doesn't support" - " a 'shape' parameter" - ) - - def expect(k: int) -> float: - return rv.dist.expect( - lambda x: x**k, args=shape, loc=loc, scale=scale - ) - - raw_moments = [expect(k) for k in range(1, 9)] - moments = raw2central(list(raw_moments)) - for k in range(2, 8): - moments[k] /= sd ** (k + 1) - - else: - if rv.dist.numargs != 0: - raise ValueError( - "The distribution provided requires a third" - " 'shape' parameter" - ) - - def expect(k: int) -> float: - return rv.dist.expect(lambda x: x**k) - - raw_moments = [expect(k) for k in range(1, 9)] - moments = raw2central(list(raw_moments)) + # Integrate the standardized central moments directly rather than + # going via raw moments and differencing them. ``loc`` and + # ``scale`` must be passed through, and the differencing route + # loses nearly all significant digits once the mean is large + # compared with the spread (e.g. N(24, 1), where the 8th raw + # moment is ~1e11 but the 8th central moment is ~1e2). + def standardized_central(k: int) -> float: + return rv.dist.expect( + lambda x: ((x - mn) / sd) ** k, + args=shape, + loc=loc, + scale=scale, + ) - moments[0] = mn # mean - moments[1] = sd**2 # variance + moments = [mn, sd**2] + [ + standardized_central(k) for k in range(3, 9) + ] self._dist = rv diff --git a/soerp_examples.py b/soerp_examples.py index 04ec828..6e0c090 100644 --- a/soerp_examples.py +++ b/soerp_examples.py @@ -20,8 +20,8 @@ print("Results should be about:") print(" > Mean................... 1176.45") print(" > Variance............... 99699.682") -print(" > Skewness Coefficient... 0.70801305") -print(" > Kurtosis Coefficient... 6.1632855") +print(" > Skewness Coefficient... 0.71007690") +print(" > Kurtosis Coefficient... 6.1611612") print(Z) print("*" * 80) @@ -87,8 +87,8 @@ print("Results should be about:") print(" > Mean................... 1176.45") print(" > Variance............... 99699.682") - print(" > Skewness Coefficient... 0.70801305") - print(" > Kurtosis Coefficient... 6.1632855") + print(" > Skewness Coefficient... 0.71007690") + print(" > Kurtosis Coefficient... 6.1611612") print(Z) print("*" * 80) @@ -184,8 +184,8 @@ print("Results should be about:") print(" > Mean................... 1176.45") print(" > Variance............... 99699.682") - print(" > Skewness Coefficient... 0.70801305") - print(" > Kurtosis Coefficient... 6.1632855") + print(" > Skewness Coefficient... 0.71007690") + print(" > Kurtosis Coefficient... 6.1611612") print(Z) print("*" * 80) diff --git a/tests/test_assembly.py b/tests/test_assembly.py index 63054f0..ff24c41 100644 --- a/tests/test_assembly.py +++ b/tests/test_assembly.py @@ -26,8 +26,13 @@ def _moments(uf): class TestThreePartAssembly: MEAN = 1176.45 VAR = 99699.682 - SKEW = 0.70801305 - KURT = 6.1632855 + # Cox (1979) Table 2 prints 81.25 / 99,699.68 / 0.7100789 / 6.18 for this + # problem. SKEW below agrees with the paper to 6 significant figures. + # KURT is the exact value of eq. (A-9) for this polynomial, cross-checked + # against an independent tensor-product Gauss quadrature; Table 2's printed + # 6.183188 is not reproducible from the paper's own equation (A-9). + SKEW = 0.71007690 + KURT = 6.1611612 def test_moments_input(self): x1 = uv([24, 1, 0, 3, 0, 15, 0, 105]) From dd7cf2d022d3812ea05b76a248db6d6a7d314bd8 Mon Sep 17 00:00:00 2001 From: saudzahirr Date: Sat, 15 Aug 2026 14:39:33 +0500 Subject: [PATCH 2/2] Fix Sonar maintainability gate and cover the new branches soerp_examples.py printed the same four "Results should be about" literals for each of the three ways every example is built, which SonarQube flagged as duplicated string literals on new code. The reference values now live in one tuple per example alongside a print_expected() helper. The script's output is byte-for-byte unchanged. The previous commit also left four lines uncovered: - umath.cot() had no test at all, so its corrected derivatives went unverified by CI. TestCot pins cot' = -csc**2 and cot'' = 2*cot*csc**2, which is precisely what the earlier tan()-based version got wrong. It sits in its own class to stay under the per-class public-method limit. - The two shape-parameter guards in UncertainVariable.__init__ were never exercised. scipy validates shape arguments when a distribution is frozen, so an inconsistent rv cannot be constructed directly; the tests set `args` on the frozen object to reach soerp's own guard, which is also what protects duck-typed, non-scipy rv objects. tests/test_uncertain_variable.py additionally locks down the standardization fix from the previous commit: moments 3-8 must depend only on distribution shape, so U(2, 6) now yields mu_4 = 1.8 and the normal moments are invariant to loc and scale. Co-authored-by: gemini-code-assist Co-authored-by: Junie Co-authored-by: Qwen-Coder Co-authored-by: Cursor Agent Co-authored-by: kilocode-bot --- .sonarcloud.properties | 11 +++- soerp_examples.py | 91 +++++++++++--------------------- tests/test_umath.py | 28 ++++++++++ tests/test_uncertain_variable.py | 72 +++++++++++++++++++++++++ 4 files changed, 141 insertions(+), 61 deletions(-) create mode 100644 tests/test_uncertain_variable.py diff --git a/.sonarcloud.properties b/.sonarcloud.properties index 9cb6219..8737852 100644 --- a/.sonarcloud.properties +++ b/.sonarcloud.properties @@ -13,4 +13,13 @@ # `--no-build`. That flag cannot be used by a project built from # source: uv refuses with "can't be installed because it is marked as # `--no-build` but has no binary distribution". -sonar.exclusions=.github/**,tests/**,docs/** +# +# soerp_examples.py is a demonstration script rather than part of the +# `soerp` package. It runs each worked example three ways - from given +# moments, from scipy.stats distributions, and from the convenience +# constructors - and prints the same reference values each time. That +# repetition is the point of the script: each block must be readable on +# its own next to the result it annotates. python:S1192 reads it as +# duplicated string literals, which fails the quality gate without +# describing any risk to users. +sonar.exclusions=.github/**,tests/**,docs/**,soerp_examples.py diff --git a/soerp_examples.py b/soerp_examples.py index 6e0c090..62941ff 100644 --- a/soerp_examples.py +++ b/soerp_examples.py @@ -9,6 +9,25 @@ from soerp import Chi2, Exp, Gamma, N, umath, uv +# Reference values for each example, kept as strings so the printed output is +# byte-for-byte what the original script emitted. Each example is run three +# times (given moments, scipy distributions, convenience constructors), so +# these live here rather than being repeated at every call site. +ASSEMBLY_EXPECTED = ("1176.45", "99699.682", "0.71007690", "6.1611612") +ORIFICE_EXPECTED = ("1330.9997", "58.210763", "0.010942207", "3.0003269") +STACKUP_EXPECTED = ("4.5", "0.75", "0.385", "3.22") +SCHEDULING_EXPECTED = ("51.7", "33.3", "0.52", "3.49") + + +def print_expected(mean: str, var: str, skew: str, kurt: str) -> None: + """Print the reference values that the result below should reproduce.""" + print("Results should be about:") + print(f" > Mean................... {mean}") + print(f" > Variance............... {var}") + print(f" > Skewness Coefficient... {skew}") + print(f" > Kurtosis Coefficient... {kurt}") + + print("UNCERTAIN DISTRIBUTION TEST FUNCTIONS USING GIVEN MOMENTS") print("*" * 80) print("Example of a three part assembly") @@ -17,11 +36,7 @@ x3 = uv([0.5, 0.25, 2, 9, 44, 265, 1854, 14833]) # exponentially distributed Z = (x1 * x2**2) / (15 * (1.5 + x3)) -print("Results should be about:") -print(" > Mean................... 1176.45") -print(" > Variance............... 99699.682") -print(" > Skewness Coefficient... 0.71007690") -print(" > Kurtosis Coefficient... 6.1611612") +print_expected(*ASSEMBLY_EXPECTED) print(Z) print("*" * 80) @@ -32,11 +47,7 @@ t = uv([165, 0.25, 0, 3, 0, 15, 0, 105]) # normally distributed C = 38.4 Q = C * umath.sqrt((520 * H * P) / (M * (t + 460))) -print("Results should be about:") -print(" > Mean................... 1330.9997") -print(" > Variance............... 58.210763") -print(" > Skewness Coefficient... 0.010942207") -print(" > Kurtosis Coefficient... 3.0003269") +print_expected(*ORIFICE_EXPECTED) print(Q) print("*" * 80) @@ -45,11 +56,7 @@ y = uv([1.5, 0.25, 2 / 3.0, 11 / 3.0, 0, 0, 0, 0]) # gamma distributed z = uv([1.5, 0.25, 2 / 3.0, 11 / 3.0, 0, 0, 0, 0]) # gamma distributed w = x + y + z -print("Results should be about:") -print(" > Mean................... 4.5") -print(" > Variance............... 0.75") -print(" > Skewness Coefficient... 0.385") -print(" > Kurtosis Coefficient... 3.22") +print_expected(*STACKUP_EXPECTED) print(w) print("*" * 80) @@ -61,11 +68,7 @@ s5 = uv([0.2, 0.04, 2, 9, 0, 0, 0, 0]) # exponental distributed s6 = uv([10, 20, 0.89, 4.2, 0, 0, 0, 0]) # chi-square distributed T = s1 + s2 + s3 + s4 + s5 + s6 -print("Results should be about:") -print(" > Mean................... 51.7") -print(" > Variance............... 33.3") -print(" > Skewness Coefficient... 0.52") -print(" > Kurtosis Coefficient... 3.49") +print_expected(*SCHEDULING_EXPECTED) print(T) ############################################################################### @@ -84,11 +87,7 @@ x2 = uv(rv=ss.norm(loc=37, scale=4)) # normally distributed x3 = uv(rv=ss.expon(scale=1 / 2.0)) # exponentially distributed Z = (x1 * x2**2) / (15 * (1.5 + x3)) - print("Results should be about:") - print(" > Mean................... 1176.45") - print(" > Variance............... 99699.682") - print(" > Skewness Coefficient... 0.71007690") - print(" > Kurtosis Coefficient... 6.1611612") + print_expected(*ASSEMBLY_EXPECTED) print(Z) print("*" * 80) @@ -99,11 +98,7 @@ t = uv(rv=ss.norm(loc=165, scale=0.5)) C = 38.4 Q = C * umath.sqrt((520 * H * P) / (M * (t + 460))) - print("Results should be about:") - print(" > Mean................... 1330.9997") - print(" > Variance............... 58.210763") - print(" > Skewness Coefficient... 0.010942207") - print(" > Kurtosis Coefficient... 3.0003269") + print_expected(*ORIFICE_EXPECTED) print(Q) print("*" * 80) @@ -119,11 +114,7 @@ y = uv(rv=ss.gamma(shape, scale=scale)) z = uv(rv=ss.gamma(shape, scale=scale)) w = x + y + z - print("Results should be about:") - print(" > Mean................... 4.5") - print(" > Variance............... 0.75") - print(" > Skewness Coefficient... 0.385") - print(" > Kurtosis Coefficient... 3.22") + print_expected(*STACKUP_EXPECTED) print(w) print("*" * 80) @@ -143,11 +134,7 @@ s5 = uv(rv=ss.expon(scale=0.2)) s6 = uv(rv=ss.chi2(10)) T = s1 + s2 + s3 + s4 + s5 + s6 - print("Results should be about:") - print(" > Mean................... 51.7") - print(" > Variance............... 33.3") - print(" > Skewness Coefficient... 0.52") - print(" > Kurtosis Coefficient... 3.49") + print_expected(*SCHEDULING_EXPECTED) print(T) print("*" * 80) @@ -181,11 +168,7 @@ x2 = N(37, 4) # normally distributed x3 = Exp(2) # exponentially distributed Z = (x1 * x2**2) / (15 * (1.5 + x3)) - print("Results should be about:") - print(" > Mean................... 1176.45") - print(" > Variance............... 99699.682") - print(" > Skewness Coefficient... 0.71007690") - print(" > Kurtosis Coefficient... 6.1611612") + print_expected(*ASSEMBLY_EXPECTED) print(Z) print("*" * 80) @@ -196,11 +179,7 @@ t = N(165, 0.5) C = 38.4 Q = C * umath.sqrt((520 * H * P) / (M * (t + 460))) - print("Results should be about:") - print(" > Mean................... 1330.9997") - print(" > Variance............... 58.210763") - print(" > Skewness Coefficient... 0.010942207") - print(" > Kurtosis Coefficient... 3.0003269") + print_expected(*ORIFICE_EXPECTED) print(Q) print("*" * 80) @@ -216,11 +195,7 @@ y = Gamma(shape, scale) z = Gamma(shape, scale) w = x + y + z - print("Results should be about:") - print(" > Mean................... 4.5") - print(" > Variance............... 0.75") - print(" > Skewness Coefficient... 0.385") - print(" > Kurtosis Coefficient... 3.22") + print_expected(*STACKUP_EXPECTED) print(w) print("*" * 80) @@ -240,11 +215,7 @@ s5 = Exp(5) s6 = Chi2(10) T = s1 + s2 + s3 + s4 + s5 + s6 - print("Results should be about:") - print(" > Mean................... 51.7") - print(" > Variance............... 33.3") - print(" > Skewness Coefficient... 0.52") - print(" > Kurtosis Coefficient... 3.49") + print_expected(*SCHEDULING_EXPECTED) print(T) print("*" * 80) diff --git a/tests/test_umath.py b/tests/test_umath.py index e982a2c..4a6d1d5 100644 --- a/tests/test_umath.py +++ b/tests/test_umath.py @@ -92,3 +92,31 @@ def test_ceil_uncertain(self): x = N(3.2, 0.01) r = umath.ceil(x) assert r.mean == pytest.approx(4.0) + + +class TestCot: + r"""cot() is separated out because its derivatives are easy to get wrong: + cot' = $-\csc^2$ and cot'' = $2\cot\csc^2$, and $\csc$ is not $\tan$.""" + + X0 = 0.7 + + def test_cot_scalar(self): + assert umath.cot(self.X0) == pytest.approx(1.0 / math.tan(self.X0)) + + def test_cot_derivatives(self): + """Regression guard: cot() used to build tan() where csc() was + required, which made both derivatives wrong.""" + x = N(self.X0, 1e-8) + r = umath.cot(x) + cot0 = math.cos(self.X0) / math.sin(self.X0) + csc2 = 1.0 / math.sin(self.X0) ** 2 + assert r.x == pytest.approx(cot0) + assert r.d(x) == pytest.approx(-csc2) + assert r.d2(x) == pytest.approx(2.0 * cot0 * csc2) + + def test_cot_uncertain(self): + x = N(self.X0, 0.01) + r = umath.cot(x) + cot0 = math.cos(self.X0) / math.sin(self.X0) + assert r.mean == pytest.approx(cot0, rel=1e-2) + assert r.var > 0 diff --git a/tests/test_uncertain_variable.py b/tests/test_uncertain_variable.py new file mode 100644 index 0000000..dbccb91 --- /dev/null +++ b/tests/test_uncertain_variable.py @@ -0,0 +1,72 @@ +""" +Tests for building UncertainVariable objects, especially from scipy.stats +distributions, where the eight standardized moments are derived by integration. +""" + +import pytest +import scipy.stats as ss + +from soerp import U, uv + + +class TestConstructionErrors: + def test_requires_moments_or_rv(self): + with pytest.raises(ValueError, match="Either the moments"): + uv() + + # scipy validates shape arguments when a distribution is frozen, so an + # inconsistent rv cannot be built directly. ``args`` is a public attribute + # of the frozen object, so we set it afterwards to reach soerp's own guard + # (which also protects against non-scipy, duck-typed rv objects). + + def test_shape_arg_rejected_for_shapeless_distribution(self): + rv = ss.norm(loc=0, scale=1) # norm takes no shape parameters + rv.args = (0.5,) + with pytest.raises(ValueError, match="doesn't support"): + uv(rv=rv) + + def test_missing_shape_arg_rejected(self): + rv = ss.gamma(3, loc=0, scale=1) # gamma requires a shape parameter + rv.args = () + with pytest.raises(ValueError, match="requires a third"): + uv(rv=rv) + + +class TestStandardizedMoments: + """Moments 3-8 are standardized, so they depend only on distribution + shape - never on location or scale.""" + + def test_uniform_moments_are_standardized(self): + """Regression: the branch for distributions without a shape parameter + ignored loc/scale and skipped standardization entirely, so U(2, 6) + reported mu_4 = 0.0125 instead of 1.8. + + For a standardized uniform, mu_(2k) = 3**k / (2k + 1). + """ + x = U(2, 6) + mn, vr, *standardized = x.moments() + assert mn == pytest.approx(4.0) + assert vr == pytest.approx(4.0 / 3.0) + assert standardized == pytest.approx( + [0.0, 9 / 5, 0.0, 27 / 7, 0.0, 81 / 9], abs=1e-6 + ) + + @pytest.mark.parametrize( + ("loc", "scale"), [(0.0, 1.0), (24.0, 1.0), (-5.0, 4.0), (1000.0, 0.5)] + ) + def test_normal_moments_invariant_to_loc_and_scale(self, loc, scale): + x = uv(rv=ss.norm(loc=loc, scale=scale)) + mn, vr, *standardized = x.moments() + assert mn == pytest.approx(loc) + assert vr == pytest.approx(scale**2) + assert standardized == pytest.approx( + [0.0, 3.0, 0.0, 15.0, 0.0, 105.0], abs=1e-6 + ) + + def test_exponential_moments(self): + x = uv(rv=ss.expon(scale=0.5)) + _, vr, *standardized = x.moments() + assert vr == pytest.approx(0.25) + assert standardized == pytest.approx( + [2.0, 9.0, 44.0, 265.0, 1854.0, 14833.0], rel=1e-6 + )