Minor updates to inverse CDF code for t skew

.gitignore

-data/

src/distribution_tools/model_imu_data.py

@@ -4,7 +4,7 @@ import scipy.stats as stats
 
 from scipy.optimize import minimize
 from tskew.tskew import getObjectiveFunction, tspdf_1d
-
+from tskew.tskew import ts_invcdf
 
 import pyreadr
 
@@ -30,16 +30,26 @@ df = 1000
 skew = 0
 
 theta = np.array([loc, scale, df, skew])
-
-res = minimize(getObjectiveFunction(realization, use_loglikelihood=True), x0=theta,
-               method='Nelder-Mead')
+#
+# res = minimize(getObjectiveFunction(realization, use_loglikelihood=True), x0=theta,
+#                method='Nelder-Mead')
 
 N = 1_000
-extent =  np.max(realization) - np.min(realization)
-xvals = np.linspace(np.min(realization) - 0.1 * extent, np.max(realization) + 0.1 * extent, N)
-
-plt.figure()
-est_pdf = tspdf_1d(xvals, res.x[0], res.x[1], res.x[2], res.x[3])
-plt.hist(realization, bins=100, density=True)
-plt.plot(xvals, est_pdf, linestyle='--', label='Estimated skew t', linewidth=4)
-plt.legend()
\ No newline at end of file
+xmin = -0.05
+xmax = 0.05
+extent =  xmax - xmin
+xvals = np.linspace(xmin - 0.1 * extent, xmax + 0.1 * extent, N)
+
+# plt.figure()
+# est_pdf = tspdf_1d(xvals, res.x[0], res.x[1], res.x[2], res.x[3])
+# plt.hist(realization, bins=500, density=True, color='green', alpha=0.5)
+# plt.plot(xvals, est_pdf, linestyle='--', label='Estimated skew t', linewidth=3, alpha=0.5)
+# plt.xlim([xmin, xmax])
+# plt.legend()
+
+loc = 2.72e-5
+scale = 2.25e-6
+df = 1
+skew = 2.8e-3
+
+# median = ts_invcdf(np.array([0.25, 0.5, 0.75]), loc, scale, df, skew)
\ No newline at end of file

src/distribution_tools/tskew/tskew.py

@@ -18,9 +18,8 @@ from scipy.optimize import minimize
 import warnings
 import scipy.io as sio
 import scipy
-# from root_finding import newton, brentq
+from .root_finding import newton, brentq
 
-from scipy.special import gamma
 
 
 plt.close('all')
@@ -28,20 +27,6 @@ from numba import vectorize, njit
 import numba as nb
 
 
-
-
-# Source: https://gregorygundersen.com/blog/2020/01/20/multivariate-t/
-# @article{kollo2021multivariate,
-#   title={Multivariate Skew t-Distribution: Asymptotics for Parameter Estimators and Extension to Skew t-Copula},
-#   author={Kollo, T{\~o}nu and K{\"a}{\"a}rik, Meelis and Selart, Anne},
-#   journal={Symmetry},
-#   volume={13},
-#   number={6},
-#   pages={1059},
-#   year={2021},
-#   publisher={Multidisciplinary Digital Publishing Institute}
-# }
-
 @njit
 def tcdf_1d(x, df):
     # if use_scipy:
@@ -234,32 +219,52 @@ def root_Newton_Rhapson(fun, x0, jac, tol=1e-12, maxiter=100):
 
     return x1
 
-def ts_invcdf(q, loc, scale, df, skew):
+def ts_invcdf_opt(q, loc, scale, df, skew):
 
     def ffun(x):
         return tscdf(x, loc, scale, df, skew) - q
 
     def fprime(x):
         return tspdf_1d(x, loc, scale, df, skew)
+
     # Do a single Newton iteration
     p0 = 0.0
     fval = ffun(p0)
     fder = fprime(p0)
     newton_step = fval / fder
+
+
     # Newton step
     p = p0 - newton_step
 
+    while ffun(p) * fval > 0:
+        p = p - newton_step
+
     if p0 < p:
-        r = brentq(f, p0, p)
+        r = brentq(ffun, p0, p)
     else:
-        r = brentq(f, p, p0)
+        r = brentq(ffun, p, p0)
 
     xvals = np.linspace(-20, 20, 1_000)
-    fvals = f(xvals)
+    fvals = ffun(xvals)
     # r = newton(func = f, x0 = 0.0, fprime = fprime)
 
     return r
 
+
+def ts_invcdf(quantiles, loc, scale, df, skew):
+    try:
+        roots = np.zeros_like(quantiles)
+        for count, q in enumerate(quantiles):
+            r = ts_invcdf_opt(q, loc, scale, df, skew)
+            pass
+
+    except:
+        pass
+
+    pass
+
+
 def numerical_inverse(rv_domain, cdf_vals):
     return scipy.interpolate.interp1d(cdf_vals, rv_domain, kind='cubic', fill_value="extrapolate")
     # def inv_fn(x):