Skip to content

Fix the second-order moment equations to match Cox (1979) - #52

Merged
saudzahirr merged 2 commits into
masterfrom
fix/cox-1979-moment-equations
Aug 15, 2026
Merged

Fix the second-order moment equations to match Cox (1979)#52
saudzahirr merged 2 commits into
masterfrom
fix/cox-1979-moment-equations

Conversation

@saudzahirr

Copy link
Copy Markdown
Collaborator

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.
  3. (A-9) triple sum used 6 b_ii^2 b_jk instead of 6 b_ii^2 b_jk^2.
  4. (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.
  5. (A-9) emitted two mu_i2 mu_j3 mu_k3 blocks, double counting three terms; the paper has one.
  6. (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.

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 <noreply@anthropic.com>
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
@saudzahirr saudzahirr self-assigned this Aug 15, 2026
@codecov-commenter

codecov-commenter commented Aug 15, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

@saudzahirr saudzahirr added the bug label Aug 15, 2026
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 <gemini-code-assist@users.noreply.github.com>
Co-authored-by: Junie <junie@jetbrains.com>
Co-authored-by: Qwen-Coder <qwen-coder@alibabacloud.com>
Co-authored-by: Cursor Agent <cursoragent@cursor.com>
Co-authored-by: kilocode-bot <github-actions@github.com>
@saudzahirr
saudzahirr force-pushed the fix/cox-1979-moment-equations branch from 2aeedd9 to dd7cf2d Compare August 15, 2026 10:14
@sonarqubecloud

Copy link
Copy Markdown

@saudzahirr

Copy link
Copy Markdown
Collaborator Author

@kilocode-bot review this!

@saudzahirr
saudzahirr merged commit a52dc5e into master Aug 15, 2026
14 checks passed
@saudzahirr
saudzahirr deleted the fix/cox-1979-moment-equations branch August 15, 2026 10:18
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants