Reorganize files.

python/# BetaPrime vs Gamma, hard-mono fit WITH.py

@@ -1,187 +0,0 @@
-# Constrained Bayesian fit (Gamma for NC, Beta-Prime for AE) — no regularization, no r-prior
-import numpy as np
-import matplotlib.pyplot as plt
-from scipy import optimize
-from scipy.special import betaln, gammaln
-import scipy.io as io
-
-data_path = "../data/"
-suv   = io.loadmat(data_path + "suv_percentilesSLOthenUWM.mat")['lung_SUVperc_COMBINED'][0:58, :, :]
-flags = io.loadmat(data_path + "flags_combined.mat")['flags'][0:58, 3]   # 0=NC, 1=AE
-
-# Feature X = max SUV_94 per subject; label y = flags
-X = np.nanmax(suv[:, :, 94], axis=1).astype(float).ravel()
-y = np.asarray(flags, int).ravel()
-
-# Guard for logs
-X = np.clip(X, 1e-12, None)
-p_emp = float(y.mean())
-
-# Small helpers
-
-def logistic(z):
-    z = np.clip(z, -60, 60)
-    return 1.0 / (1.0 + np.exp(-z))
-
-def sigmoid(t):
-    return 1.0 / (1.0 + np.exp(-t))
-
-def softplus(t):
-    t = np.asarray(t, float)
-    return np.log1p(np.exp(-np.abs(t))) + np.maximum(t, 0.0)
-
-# dE(x) pieces for log-odds
-def dE_ess(x, a, b, s, k, th):
-    x = np.asarray(x, float)
-    return (a - k) * np.log(x) - (a + b) * np.log1p(x / s) + x / th
-
-def dE_const(a, b, s, k, th):
-    return -(a * np.log(s)) - betaln(a, b) + k * np.log(th) + gammaln(k)
-
-def dE_full(x, a, b, s, k, th):
-    return dE_ess(x, a, b, s, k, th) + dE_const(a, b, s, k, th)
-
-# Monotonicity cap for theta
-def theta_max(a, b, k, s, eps=1e-12):
-    A = a - k
-    if A <= 0:
-        return np.inf
-    r = np.sqrt(a + b) - np.sqrt(max(A, eps))
-    if r <= 1e-12:
-        return np.inf
-    return s / (r * r)
-
-# φ = [p_raw, b_raw, s_raw, k_raw, d_raw, u_raw]  (unconstrained)
-def unpack_phi_mono(phi):
-    p_raw, b_raw, s_raw, k_raw, d_raw, u_raw = phi
-    p = sigmoid(p_raw)                       # (0,1)
-    b = softplus(b_raw) + 1e-6               # >0
-    s = softplus(s_raw) + 1e-6               # >0
-    k = softplus(k_raw) + 1e-6               # >0
-    delta = softplus(d_raw) + 1e-6           # >0
-    a = k + delta                            # enforce a > k
-    th_cap = theta_max(a, b, k, s)           # theta cap from monotonicity
-    th = th_cap * sigmoid(u_raw)             # 0 < theta <= th_cap
-    return p, a, b, s, k, th
-
-#  Prior on p 
-TAU = 25.0                                   # shrink toward empirical AE rate
-alpha = max(TAU * p_emp, 1e-6)
-beta  = max(TAU * (1.0 - p_emp), 1e-6)
-
-prior_r = None
-#prior_r = (1.01, 1.01)
-#prior_r = (3, 3)
-
-# Objective: negative log-posterior (likelihood + Beta prior on p)
-def neg_post_phi_mono(phi, X, y):
-    p, a, b, s, k, th = unpack_phi_mono(phi)
-    eps = 1e-12
-
-    L  = (np.log(p) - np.log(1 - p)) + dE_full(X, a, b, s, k, th)
-    px = logistic(L)
-    nll = -np.sum(y * np.log(px + eps) + (1 - y) * np.log(1 - px + eps))
-
-    # Beta(alpha, beta) prior on p → negative log-prior
-    npr_p = -((alpha - 1) * np.log(p + eps) + (beta - 1) * np.log(1 - p + eps))
-
-    if prior_r is None: return nll + npr_p
-        
-    # Prior on r = theta / theta_max (softly avoid boundaries)
-    thcap = theta_max(a, b, k, s)
-    if np.isfinite(thcap) and thcap > 0:
-        r = np.clip(th / thcap, 1e-9, 1 - 1e-9)
-        # Negative log Beta prior:  -[(α-1)log r + (β-1)log(1-r)]  (const dropped)
-        npr_r = -((prior_r[0] -1) * np.log(r) + (prior_r[1]-1)* np.log(1.0 - r))
-    else:
-        npr_r = 0.0
-
-    return nll + npr_p + npr_r
-
-# Initialization (stable, simple)
-def init_phi(X, y):
-    # Gamma(k, theta) MoM for NC group (only need k0 as a safe size proxy)
-    X0 = X[y == 0]
-    m0 = X0.mean() if X0.size else X.mean()
-    v0 = X0.var()  if X0.size else X.var()
-    k0 = 2.0 if v0 <= 0 else max((m0**2)/(v0 + 1e-9), 1.5)
-
-    # AE median to seed s0
-    X1 = X[y == 1]
-    m1 = np.median(X1) if X1.size else np.median(X)
-
-    p0 = np.clip(float(y.mean()), 1e-3, 1 - 1e-3)
-    b0, s0 = 1.5, max(m1, 0.5)
-
-    return np.array([
-        np.log(p0 / (1 - p0)),             # p_raw
-        np.log(np.expm1(b0) + 1e-9),       # b_raw
-        np.log(np.expm1(s0) + 1e-9),       # s_raw
-        np.log(np.expm1(k0) + 1e-9),       # k_raw
-        np.log(np.expm1(1.0) + 1e-9),      # d_raw  (delta)
-        -0.2                               # u_raw  (keeps theta a bit below cap initially)
-    ], float)
-
-#  Fit wrapper (one retry)
-def fit_bayes_mono(X, y, phi_start=None, rng=None):
-    if rng is None:
-        rng = np.random.default_rng(0)
-    if phi_start is None:
-        phi_start = init_phi(X, y)
-
-    obj = lambda phi: neg_post_phi_mono(phi, X, y)
-
-    res = optimize.minimize(
-        obj, phi_start, method="L-BFGS-B",
-        options={"maxiter": 6000, "ftol": 1e-9}
-    )
-    if not (res.success and np.isfinite(res.fun)):
-        phi_try = phi_start + rng.normal(0, 0.2, size=phi_start.shape)
-        res = optimize.minimize(
-            obj, phi_try, method="L-BFGS-B",
-            options={"maxiter": 6000, "ftol": 1e-9}
-        )
-    return unpack_phi_mono(res.x), res
-
-#  prediction 
-def P_with(theta, x):
-    p, a, b, s, k, th = theta
-    L = (np.log(p) - np.log(1 - p)) + dE_full(x, a, b, s, k, th)
-    return logistic(L)
-
-# Run fit + plot
-theta_hat, res = fit_bayes_mono(X, y)
-print("Optimization success:", res.success, " fval:", float(res.fun))
-
-(p,a,b,s,k,th) = theta_hat
-thcap = theta_max(a, b, k, s)
-print("theta (p,a,b,s,k,theta):", tuple(float(t) for t in theta_hat), "ratio(th):", th/thcap)
-
-# x-range (cap right end at 10 for readability)
-x_lo = max(1e-6, float(X.min()) * 0.8)
-x_hi = min(10.0, float(X.max()) * 1.2)
-xg   = np.linspace(x_lo, x_hi, 600)
-p_curve = P_with(theta_hat, xg)
-
-fig, ax = plt.subplots(figsize=(7.0, 4.6), dpi=140)
-ax.plot(xg, p_curve, color="#000000", lw=2.2, label="P(AE|x) (MAP)")
-
-# overlay data with tiny vertical jitter so points don't overlap
-rng_plot = np.random.default_rng(999)
-jit = (rng_plot.random(len(y)) - 0.5) * 0.06
-ax.scatter(X[y==0], (y + jit)[y==0], s=22, alpha=0.55, color="#2ca02c", edgecolors='none', label='NC')
-ax.scatter(X[y==1], (y + jit)[y==1], s=26, alpha=0.75, color="#ff7f0e", edgecolors='none', label='AE')
-
-ax.set_ylim(-0.05, 1.05)
-ax.set_xlabel('x')
-ax.set_ylabel('P(AE | x)')
-
-if prior_r is None:
-    ax.set_title('Constrained Bayesian fit (no regularization, prior on p)')
-else:
-    ax.set_title(f'Constrained Bayesian fit (no regularization, prior on p and prior r{prior_r})')
-
-ax.grid(alpha=0.3)
-ax.legend(loc='lower right', frameon=False)
-plt.tight_layout()
-plt.show()

python/Untitled-1.ipynb

@@ -1,208 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "eae6df92",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Constrained Bayesian fit (Gamma for NC, Beta-Prime for AE) — no regularization, no r-prior\n",
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt\n",
-    "from scipy import optimize\n",
-    "from scipy.special import betaln, gammaln\n",
-    "import scipy.io as io\n",
-    "\n",
-    "data_path = \"../data/\"\n",
-    "suv   = io.loadmat(data_path + \"suv_percentilesSLOthenUWM.mat\")['lung_SUVperc_COMBINED'][0:58, :, :]\n",
-    "flags = io.loadmat(data_path + \"flags_combined.mat\")['flags'][0:58, 3]   # 0=NC, 1=AE\n",
-    "\n",
-    "# Feature X = max SUV_94 per subject; label y = flags\n",
-    "X = np.nanmax(suv[:, :, 94], axis=1).astype(float).ravel()\n",
-    "y = np.asarray(flags, int).ravel()\n",
-    "\n",
-    "# Guard for logs\n",
-    "X = np.clip(X, 1e-12, None)\n",
-    "p_emp = float(y.mean())\n",
-    "\n",
-    "# Small helpers\n",
-    "\n",
-    "def logistic(z):\n",
-    "    z = np.clip(z, -60, 60)\n",
-    "    return 1.0 / (1.0 + np.exp(-z))\n",
-    "\n",
-    "def sigmoid(t):\n",
-    "    return 1.0 / (1.0 + np.exp(-t))\n",
-    "\n",
-    "def softplus(t):\n",
-    "    t = np.asarray(t, float)\n",
-    "    return np.log1p(np.exp(-np.abs(t))) + np.maximum(t, 0.0)\n",
-    "\n",
-    "# dE(x) pieces for log-odds\n",
-    "def dE_ess(x, a, b, s, k, th):\n",
-    "    x = np.asarray(x, float)\n",
-    "    return (a - k) * np.log(x) - (a + b) * np.log1p(x / s) + x / th\n",
-    "\n",
-    "def dE_const(a, b, s, k, th):\n",
-    "    return -(a * np.log(s)) - betaln(a, b) + k * np.log(th) + gammaln(k)\n",
-    "\n",
-    "def dE_full(x, a, b, s, k, th):\n",
-    "    return dE_ess(x, a, b, s, k, th) + dE_const(a, b, s, k, th)\n",
-    "\n",
-    "# Monotonicity cap for theta\n",
-    "def theta_max(a, b, k, s, eps=1e-12):\n",
-    "    A = a - k\n",
-    "    if A <= 0:\n",
-    "        return np.inf\n",
-    "    r = np.sqrt(a + b) - np.sqrt(max(A, eps))\n",
-    "    if r <= 1e-12:\n",
-    "        return np.inf\n",
-    "    return s / (r * r)\n",
-    "\n",
-    "# φ = [p_raw, b_raw, s_raw, k_raw, d_raw, u_raw]  (unconstrained)\n",
-    "def unpack_phi_mono(phi):\n",
-    "    p_raw, b_raw, s_raw, k_raw, d_raw, u_raw = phi\n",
-    "    p = sigmoid(p_raw)                       # (0,1)\n",
-    "    b = softplus(b_raw) + 1e-6               # >0\n",
-    "    s = softplus(s_raw) + 1e-6               # >0\n",
-    "    k = softplus(k_raw) + 1e-6               # >0\n",
-    "    delta = softplus(d_raw) + 1e-6           # >0\n",
-    "    a = k + delta                            # enforce a > k\n",
-    "    th_cap = theta_max(a, b, k, s)           # theta cap from monotonicity\n",
-    "    th = th_cap * sigmoid(u_raw)             # 0 < theta <= th_cap\n",
-    "    return p, a, b, s, k, th\n",
-    "\n",
-    "#  Prior on p \n",
-    "TAU = 25.0                                   # shrink toward empirical AE rate\n",
-    "alpha = max(TAU * p_emp, 1e-6)\n",
-    "beta  = max(TAU * (1.0 - p_emp), 1e-6)\n",
-    "\n",
-    "prior_r = None\n",
-    "#prior_r = (1.01, 1.01)\n",
-    "#prior_r = (1.05, 1.05)\n",
-    "#prior_r = (3, 3)\n",
-    "#prior_r = (1.2, 1.2)\n",
-    "# Objective: negative log-posterior (likelihood + Beta prior on p)\n",
-    "def neg_post_phi_mono(phi, X, y):\n",
-    "    p, a, b, s, k, th = unpack_phi_mono(phi)\n",
-    "    eps = 1e-12\n",
-    "\n",
-    "    L  = (np.log(p) - np.log(1 - p)) + dE_full(X, a, b, s, k, th)\n",
-    "    px = logistic(L)\n",
-    "    nll = -np.sum(y * np.log(px + eps) + (1 - y) * np.log(1 - px + eps))\n",
-    "\n",
-    "    # Beta(alpha, beta) prior on p → negative log-prior\n",
-    "    npr_p = -((alpha - 1) * np.log(p + eps) + (beta - 1) * np.log(1 - p + eps))\n",
-    "\n",
-    "    if prior_r is None: return nll + npr_p\n",
-    "        \n",
-    "    # Prior on r = theta / theta_max (softly avoid boundaries)\n",
-    "    thcap = theta_max(a, b, k, s)\n",
-    "    if np.isfinite(thcap) and thcap > 0:\n",
-    "        r = np.clip(th / thcap, 1e-9, 1 - 1e-9)\n",
-    "        # Negative log Beta prior:  -[(α-1)log r + (β-1)log(1-r)]  (const dropped)\n",
-    "        npr_r = -((prior_r[0] -1) * np.log(r) + (prior_r[1]-1)* np.log(1.0 - r))\n",
-    "    else:\n",
-    "        npr_r = 0.0\n",
-    "\n",
-    "    return nll + npr_p + npr_r\n",
-    "\n",
-    "# Initialization (stable, simple)\n",
-    "def init_phi(X, y):\n",
-    "    # Gamma(k, theta) MoM for NC group (only need k0 as a safe size proxy)\n",
-    "    X0 = X[y == 0]\n",
-    "    m0 = X0.mean() if X0.size else X.mean()\n",
-    "    v0 = X0.var()  if X0.size else X.var()\n",
-    "    k0 = 2.0 if v0 <= 0 else max((m0**2)/(v0 + 1e-9), 1.5)\n",
-    "\n",
-    "    # AE median to seed s0\n",
-    "    X1 = X[y == 1]\n",
-    "    m1 = np.median(X1) if X1.size else np.median(X)\n",
-    "\n",
-    "    p0 = np.clip(float(y.mean()), 1e-3, 1 - 1e-3)\n",
-    "    b0, s0 = 1.5, max(m1, 0.5)\n",
-    "\n",
-    "    return np.array([\n",
-    "        np.log(p0 / (1 - p0)),             # p_raw\n",
-    "        np.log(np.expm1(b0) + 1e-9),       # b_raw\n",
-    "        np.log(np.expm1(s0) + 1e-9),       # s_raw\n",
-    "        np.log(np.expm1(k0) + 1e-9),       # k_raw\n",
-    "        np.log(np.expm1(1.0) + 1e-9),      # d_raw  (delta)\n",
-    "        -0.2                               # u_raw  (keeps theta a bit below cap initially)\n",
-    "    ], float)\n",
-    "\n",
-    "#  Fit wrapper (one retry)\n",
-    "def fit_bayes_mono(X, y, phi_start=None, rng=None):\n",
-    "    if rng is None:\n",
-    "        rng = np.random.default_rng(0)\n",
-    "    if phi_start is None:\n",
-    "        phi_start = init_phi(X, y)\n",
-    "\n",
-    "    obj = lambda phi: neg_post_phi_mono(phi, X, y)\n",
-    "\n",
-    "    res = optimize.minimize(\n",
-    "        obj, phi_start, method=\"L-BFGS-B\",\n",
-    "        options={\"maxiter\": 6000, \"ftol\": 1e-9}\n",
-    "    )\n",
-    "    if not (res.success and np.isfinite(res.fun)):\n",
-    "        phi_try = phi_start + rng.normal(0, 0.2, size=phi_start.shape)\n",
-    "        res = optimize.minimize(\n",
-    "            obj, phi_try, method=\"L-BFGS-B\",\n",
-    "            options={\"maxiter\": 6000, \"ftol\": 1e-9}\n",
-    "        )\n",
-    "    return unpack_phi_mono(res.x), res\n",
-    "\n",
-    "#  prediction \n",
-    "def P_with(theta, x):\n",
-    "    p, a, b, s, k, th = theta\n",
-    "    L = (np.log(p) - np.log(1 - p)) + dE_full(x, a, b, s, k, th)\n",
-    "    return logistic(L)\n",
-    "\n",
-    "# Run fit + plot\n",
-    "theta_hat, res = fit_bayes_mono(X, y)\n",
-    "print(\"Optimization success:\", res.success, \" fval:\", float(res.fun))\n",
-    "\n",
-    "(p,a,b,s,k,th) = theta_hat\n",
-    "thcap = theta_max(a, b, k, s)\n",
-    "print(\"theta (p,a,b,s,k,theta):\", tuple(float(t) for t in theta_hat), \"ratio(th):\", th/thcap)\n",
-    "\n",
-    "# x-range (cap right end at 10 for readability)\n",
-    "x_lo = max(1e-6, float(X.min()) * 0.8)\n",
-    "x_hi = min(10.0, float(X.max()) * 1.2)\n",
-    "xg   = np.linspace(x_lo, x_hi, 600)\n",
-    "p_curve = P_with(theta_hat, xg)\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(7.0, 4.6), dpi=140)\n",
-    "ax.plot(xg, p_curve, color=\"#000000\", lw=2.2, label=\"P(AE|x) (MAP)\")\n",
-    "\n",
-    "# overlay data with tiny vertical jitter so points don't overlap\n",
-    "rng_plot = np.random.default_rng(999)\n",
-    "jit = (rng_plot.random(len(y)) - 0.5) * 0.06\n",
-    "ax.scatter(X[y==0], (y + jit)[y==0], s=22, alpha=0.55, color=\"#2ca02c\", edgecolors='none', label='NC')\n",
-    "ax.scatter(X[y==1], (y + jit)[y==1], s=26, alpha=0.75, color=\"#ff7f0e\", edgecolors='none', label='AE')\n",
-    "\n",
-    "ax.set_ylim(-0.05, 1.05)\n",
-    "ax.set_xlabel('x')\n",
-    "ax.set_ylabel('P(AE | x)')\n",
-    "\n",
-    "if prior_r is None:\n",
-    "    ax.set_title('Constrained Bayesian fit (no regularization, prior on p)')\n",
-    "else:\n",
-    "    ax.set_title(f'Constrained Bayesian fit (no regularization, prior on p and prior r{prior_r})')\n",
-    "\n",
-    "ax.grid(alpha=0.3)\n",
-    "ax.legend(loc='lower right', frameon=False)\n",
-    "plt.tight_layout()\n",
-    "plt.show()\n"
-   ]
-  }
- ],
- "metadata": {
-  "language_info": {
-   "name": "python"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}

python/logistic/data_utils.py

@@ -1,3 +1,9 @@
+"""
+    Utility functions for data extraction and processing
+
+    Authors: Martin Horvat, January 2026
+"""
+
 import numpy as np
 
 """

python/logistic/logit_reg_fit.ipynb

@@ -1221,11 +1221,7 @@
  ],
  "metadata": {
   "kernelspec": {
-<<<<<<< HEAD
-   "display_name": "pymc-env (3.10.11)",
-=======
    "display_name": "base (3.12.3)",
->>>>>>> 4f352e73131a6474f93318ebafa21031df5e2fd0
    "language": "python",
    "name": "python3"
   },
@@ -1239,7 +1235,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.11"
+   "version": "3.12.3"
   }
  },
  "nbformat": 4,

python/logistic/logit_utils.py

@@ -1,3 +1,19 @@
+"""
+    Logistic regression utilities with polynomial logit (log odds) function:
+
+        p(x|b) = 1/(1 + exp(-F(x|b)))
+    
+    with log odds F of polynomial form:
+
+        F(x|b) = sum_{i=0}^degree b_i x^i
+                
+    with decision function (aka logit) F and parameters 
+        
+        b = [b_i]_{i=0}^degree
+
+    Author: Martin Horvat, January 2026
+"""
+
 import numpy as np
 import scipy
 import scipy.stats

python/logistic/logit_utils_gen.py

@@ -1,3 +1,29 @@
+"""
+    Providing class for logistic regression utilities with polynomial 
+    logit (log odds) function:
+
+        p(x|pars) = 1/(1 + exp(-F(x|b(pars))))
+    
+    with log odds F of polynomial form:
+
+        F(x|b) = sum_{i=0}^degree b_i x^i
+                
+    with decision function (aka logit) F and parameters 
+        
+        b(pars) = [b(pars)_i]_{i=0}^degree
+    
+    where pars are regression parameters. Coefficients
+    b(pars) can be constrained to be monotonic function of x by using 
+    monotonic cubic transformation:
+        
+        b(pars) = mc.forward_map(pars)
+    
+    where mc is module mono_cubic2.
+
+    Author: Martin Horvat, January 2026
+"""
+
+
 import numpy as np
 import scipy
 import scipy.optimize

python/logistic/mvn.py

@@ -1,15 +1,16 @@
+"""
+    Multivariate normality tests
+    
+    Authors: Martin Horvat, Janury 2026
+"""
+
 import numpy as np
 import scipy
 
 from typing import Tuple
 
-
-# * * * STATISTICAL TESTS * * *
-
 def mardia_test(data: np.ndarray, cov: bool = True) -> Tuple[float, float, float, float]:
     """
-    https://rdrr.io/cran/MVN/src/R/mvn.R
-    https://stats.stackexchange.com/questions/317147/how-to-get-a-single-p-value-from-the-two-p-values-of-a-mardias-multinormality-t
     Mardia's multivariate skewness and kurtosis.
     Calculates the Mardia's multivariate skewness and kurtosis coefficients
     as well as their corresponding statistical test. For large sample size
@@ -18,17 +19,22 @@ def mardia_test(data: np.ndarray, cov: bool = True) -> Tuple[float, float, float
     both uncorrected and corrected skewness statistic are presented. Likewise,
     the multivariate kurtosis it is distributed as a unit-normal.
 
-     Syntax: function [Mskekur] = Mskekur(X,c,alpha)
-
-     Inputs:
-          X - multivariate data matrix [Size of matrix must be n(data)-by-p(variables)].
-          cov - boolean to whether to normalize the covariance matrix by n (c=1[default]) or by n-1 (c~=1)
-
-     Outputs:
-          - skewness test statistic
-          - kurtosis test statistic
-          - significance value for skewness
-          - significance value for kurtosis
+        Syntax: function [Mskekur] = Mskekur(X,c,alpha)
+
+    Ref:
+      * https://rdrr.io/cran/MVN/src/R/mvn.R
+      * https://stats.stackexchange.com/questions/317147/how-to-get-a-single-p-value-from-the-two-p-values-of-a-mardias-multinormality-t
+   
+    Inputs:
+        X - multivariate data matrix [Size of matrix must be n(data)-by-p(variables)].
+        cov - boolean to whether to normalize the covariance matrix by n (c=1[default]) or by n-1 (c~=1)
+
+    Outputs:
+        tuple containing:
+            skewness test statistic,
+            kurtosis test statistic,
+            significance value for skewness,
+            significance value for kurtosis
     """
     n, p = data.shape
 
@@ -82,13 +88,18 @@ def hz_test(data: np.ndarray, cov: bool = True) -> Tuple[float, float]:
     """
     Henze-Zirkler method for goodness of fit of data to a multivariate normal distribution.
     Researchers tend to use this MVN test for larger samples (N > 100).
-    https://www.tandfonline.com/doi/abs/10.1080/03610929008830400
 
-    :param data: multivariate data matrix [Size of matrix must be n(data)-by-p(variables)].
-    :param cov: boolean to whether to normalize the covariance matrix by n (c=1[default]) or by n-1 (c~=1)
-    :return:
-        HZ - Henze-Zirkler test statistic
-        p_value - significance value
+    Ref:
+    * https://www.tandfonline.com/doi/abs/10.1080/03610929008830400
+
+    Input:
+        data: multivariate data matrix [Size of matrix must be n(data)-by-p(variables)].
+        cov: boolean to whether to normalize the covariance matrix by n (c=1[default]) or by n-1 (c~=1)
+    
+    Return:
+        tuple containing:
+            HZ - Henze-Zirkler test statistic
+            p_value - significance value
     """
     n, p = data.shape
 
@@ -145,7 +156,7 @@ def royston_test(X):
     """
     Royston's Multivariate Normality Test using Fisher's method on Shapiro-Wilk p-values.
     
-    Parameters:
+    Input:
         X (ndarray): 2D array (n_samples x n_variables)
 
     Returns: