31.2 xi in #17: error band estimator

pymule/xicut.py

 from scipy import optimize
+import scipy.stats
 import numpy as np
 
 
@@ -18,3 +19,25 @@ def myfit(data, n):
     covar = out[1]
 
     return coeff, covar
+
+
+def get_val(xs, coeff):
+    return [
+        sum(np.log(x)**np.arange(0,len(coeff)) * coeff)
+        for x in xs
+    ]
+
+
+def get_errorbands(x, coeff, covar, cf=0.9):
+    if type(x) == np.ndarray:
+        return np.array([
+            get_errorbands(i, coeff, covar, cf) for i in x
+        ])
+    else:
+        tval = scipy.stats.t.ppf((cf+1)/2, len(coeff))
+        xvec = np.log(x)**np.arange(0,len(coeff))
+
+        delta = tval*np.sqrt(np.matmul(np.matmul(xvec, covar), xvec))
+        centre = sum(xvec * coeff)
+
+        return np.array([x, centre, centre+delta, centre-delta])