Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
saudzahirr merged 2 commits into from
Aug 15, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
cb2e3a5
Fix the second-order moment equations to match Cox (1979)
saudzahirr Aug 15, 2026
dd7cf2d
Fix Sonar maintainability gate and cover the new branches
saudzahirr Aug 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 10 additions & 1 deletion .sonarcloud.properties
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
8 changes: 4 additions & 4 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -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
Expand Down
60 changes: 47 additions & 13 deletions docs/theory.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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$):
Expand All @@ -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\}
$$

Expand Down
27 changes: 8 additions & 19 deletions soerp/method_of_moments.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
)

############################
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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]
Expand All @@ -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]
Expand All @@ -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]
Expand Down Expand Up @@ -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]
Expand All @@ -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]
Expand Down Expand Up @@ -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]
Expand Down
4 changes: 2 additions & 2 deletions soerp/umath.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
61 changes: 28 additions & 33 deletions soerp/uncertain_variable.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
import numpy as np
import scipy.stats as ss

from .method_of_moments import raw2central
from .uncertain_function import UncertainFunction


Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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

Expand Down
Loading
Loading