Document monotonic parameter sign symmetries

src/irae_risk/monotonic.py

@@ -5,16 +5,27 @@
 
     We define the polynomial:
           poly(x) = sum_i beta_i * x^i
-    where the coefficients beta_i are parameterized by the vector 'pars',
+    where the coefficients beta_i are parameterized by the vector 'theta',
     meaning:
-            beta_i is a function pars
+            beta_i is a function theta
+
+    IMPORTANT -- theta is not unique.  The forward map is invariant under
+
+        (C, epsilon, k1, k2) -> (C, -epsilon, k1, k2)
+        (C, epsilon, k1, k2) -> (C, epsilon, -k1, -k2).
+
+    Thus, up to four theta vectors represent the same polynomial (fewer when
+    epsilon or k1 is zero).  Likelihoods, beta coefficients, derivatives of the
+    polynomial, and fitted probabilities are invariant under these sign
+    changes.  Raw theta summaries and confidence intervals are not invariant
+    unless a canonical branch is imposed, for example epsilon >= 0 and k1 >= 0.
 """
 import numpy as np
 
 # Forward map: Parameters to polynomial coefficients
-def forward_map(pars):
+def forward_map(theta):
     """
-    Converts parameters (pars = [C, epsilon, k1, k2]) into polynomial coefficients (beta = [d, c, b, a]).
+    Converts parameters (theta = [C, epsilon, k1, k2]) into polynomial coefficients (beta = [d, c, b, a]).
     
     Polynomial Definition:
     - d = C: Constant term of the polynomial.
@@ -23,12 +34,17 @@ def forward_map(pars):
     - a = k1**2 / 3: Coefficient of the cubic term.
 
     Parameters:
-    - pars: NumPy array of parameters [C, epsilon, k1, k2].
+    - theta: NumPy array of parameters [C, epsilon, k1, k2].
 
     Returns:
     - beta: NumPy array of polynomial coefficients [d, c, b, a].
+
+    Notes:
+    - The map is many-to-one because changing the sign of epsilon, or changing
+      both signs of k1 and k2, leaves beta unchanged.  Do not interpret raw
+      theta summaries across fits without first selecting a common branch.
     """
-    C, epsilon, k1, k2 = pars  # Unpack the parameter vector
+    C, epsilon, k1, k2 = theta  # Unpack the parameter vector
     
     # Compute coefficients
     d = C
@@ -42,17 +58,17 @@ def forward_map(pars):
 
 
 # Jacobian of the forward map: First-order derivatives
-def forward_map_jacobian(pars):
+def forward_map_jacobian(theta):
     """
     Computes the Jacobian matrix of the forward map analytically.
     
     Parameters:
-    - pars: NumPy array of parameters [C, epsilon, k1, k2].
+    - theta: NumPy array of parameters [C, epsilon, k1, k2].
     
     Returns:
-    - J: NumPy 4x4 Jacobian matrix, where J[i, j] = d(beta[i])/d(pars[j]).
+    - J: NumPy 4x4 Jacobian matrix, where J[i, j] = d(beta[i])/d(theta[j]).
     """
-    C, epsilon, k1, k2 = pars  # Unpack parameters
+    C, epsilon, k1, k2 = theta  # Unpack parameters
     
     # Initialize Jacobian matrix
     J = np.zeros((4, 4))  # 4x4 matrix
@@ -85,17 +101,17 @@ def forward_map_jacobian(pars):
 
 
 # Hessian of the forward map: Second-order derivatives
-def forward_map_hessian(pars):
+def forward_map_hessian(theta):
     """
     Computes the Hessian tensor of the forward map analytically.
     
     Parameters:
-    - pars: NumPy array of parameters [C, epsilon, k1, k2].
+    - theta: NumPy array of parameters [C, epsilon, k1, k2].
     
     Returns:
-    - H: NumPy 4x4x4 Hessian tensor, where H[i, j, k] = d^2(beta[i])/d(pars[j])d(pars[k]).
+    - H: NumPy 4x4x4 Hessian tensor, where H[i, j, k] = d^2(beta[i])/d(theta[j])d(theta[k]).
     """
-    C, epsilon, k1, k2 = pars  # Unpack parameters
+    C, epsilon, k1, k2 = theta  # Unpack parameters
     
     # Initialize Hessian tensor (4 x 4 x 4)
     H = np.zeros((4, 4, 4))
@@ -120,7 +136,7 @@ def forward_map_hessian(pars):
 # Backward map: Polynomial coefficients to parameters
 def backward_map(beta, only_one = True):
     """
-    Computes the parameters (pars = [C, epsilon, k1, k2]) from the polynomial coefficients (beta = [d, c, b, a]).
+    Computes the parameters (theta = [C, epsilon, k1, k2]) from the polynomial coefficients (beta = [d, c, b, a]).
     
     Polynomial Definition:
     - d = C: Constant term of the polynomial.
@@ -174,30 +190,29 @@ def backward_map(beta, only_one = True):
 # Round-trip test
 # -------------------------
 if __name__ == "__main__":
-    pars_original = np.array([1.0, 2.0, 3.0, 4.0])
-    beta = forward_map(pars_original)
-    pars_recovered = backward_map(beta)
+    theta_original = np.array([1.0, 2.0, 3.0, 4.0])
+    beta = forward_map(theta_original)
+    theta_recovered = backward_map(beta)
 
-    print("Original pars: ", pars_original)
+    print("Original theta: ", theta_original)
     print("Beta:          ", beta)
-    print("Recovered pars:", pars_recovered)
-    print("Difference:    ", pars_recovered - pars_original)
+    print("Recovered theta:", theta_recovered)
+    print("Difference:    ", theta_recovered - theta_original)
 
-    print("\nJacobian at pars:")
-    print(forward_map_jacobian(pars_original))
+    print("\nJacobian at theta:")
+    print(forward_map_jacobian(theta_original))
 
     print("\nHessian for beta1:")
-    print(forward_map_hessian(pars_original)[1])
+    print(forward_map_hessian(theta_original)[1])
 
     print("\nLinear func:")
     lin_fun_beta = [1,0.2,0,0]
     only_one = False
 
-    lin_fun_pars = backward_map(lin_fun_beta, only_one=only_one)
-    lin_fun_beta_recover = lin_fun_pars if only_one else np.unique([forward_map(pars) for pars in lin_fun_pars], axis=0)
+    lin_fun_theta = backward_map(lin_fun_beta, only_one=only_one)
+    lin_fun_beta_recover = lin_fun_theta if only_one else np.unique([forward_map(theta) for theta in lin_fun_theta], axis=0)
     
     print(f"  {only_one = }")
     print("  lin_fun_beta:", lin_fun_beta)
-    print("  backwards:", lin_fun_pars)
+    print("  backwards:", lin_fun_theta)
     print("  forwards:", lin_fun_beta_recover)
-

tests/test_monotonic.py

@@ -0,0 +1,17 @@
+import numpy as np
+
+from irae_risk.monotonic import forward_map
+
+
+def test_forward_map_is_invariant_under_parameter_sign_symmetries():
+    theta = np.array([1.2, 0.7, -1.5, 2.3])
+    expected = forward_map(theta)
+
+    equivalent_theta = (
+        theta * np.array([1, -1, 1, 1]),
+        theta * np.array([1, 1, -1, -1]),
+        theta * np.array([1, -1, -1, -1]),
+    )
+
+    for equivalent in equivalent_theta:
+        np.testing.assert_allclose(forward_map(equivalent), expected)