Added doctest to check_random_state and added ci_quantile and ci_sf to FitDistribution

wafo/misc.py

@@ -21,7 +21,6 @@ from numdifftools.extrapolation import dea3  # @UnusedImport
 from wafo.plotbackend import plotbackend
 from collections import Callable
 import numbers
-
 try:
     from wafo import c_library as clib  # @UnresolvedImport
 except ImportError:
@@ -54,6 +53,17 @@ def check_random_state(seed):
     If seed is an int, return a new RandomState instance seeded with seed.
     If seed is already a RandomState instance, return it.
     Otherwise raise ValueError.
+
+    Example
+    -------
+    >>> check_random_state(seed=None)
+    <mtrand.RandomState object at ...
+    >>> check_random_state(seed=1)
+    <mtrand.RandomState object at ...
+    >>> check_random_state(seed=np.random.RandomState(1))
+    <mtrand.RandomState object at ...
+
+    check_random_state(seed=2.5)
     """
     if seed is None or seed is np.random:
         return np.random.mtrand._rand
@@ -61,8 +71,8 @@ def check_random_state(seed):
         return np.random.RandomState(seed)
     if isinstance(seed, np.random.RandomState):
         return seed
-    raise ValueError('%r cannot be used to seed a numpy.random.RandomState'
+    msg = '{} cannot be used to seed a numpy.random.RandomState instance'
-                     ' instance' % seed)
+    raise ValueError(msg.format(seed))
 
 
 def valarray(shape, value=np.NaN, typecode=None):
@@ -1373,8 +1383,8 @@ def mctp2tc(f_Mm, utc, param, f_mM=None):
                     if dim_p == 1:
                         tempp[0] = (Ap/(1-Bp*Ap)*e)
                     else:
-                        tempp = np.dot(Ap,
+                        rh = I - np.dot(Bp, Ap)
-                                       linalg.lstsq(I-np.dot(Bp, Ap), e)[0])
+                        tempp = np.dot(Ap, linalg.solve(rh, e))
 
                     # end
                 # end
@@ -3099,7 +3109,7 @@ def test_docstrings():
     # np.set_printoptions(precision=6)
     import doctest
     print('Testing docstrings in %s' % __file__)
-    doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
+    doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE|doctest.ELLIPSIS)
 
 
 if __name__ == "__main__":

wafo/stats/estimation.py

@@ -1345,6 +1345,26 @@ class FitDistribution(rv_frozen):
         '''
         return ProfileProbability(self, log_sf, **kwds)
 
+    def ci_sf(self, sf, alpha=0.05, i=2):
+        ci = []
+        for log_sfi in np.atleast_1d(np.log(sf)).ravel():
+            try:
+                Lp = self.profile_probability(log_sfi, i=i)
+                ci.append(np.exp(Lp.get_bounds(alpha=alpha)))
+            except Exception:
+                ci.append((np.nan, np.nan))
+        return np.array(ci)
+
+    def ci_quantile(self, x, alpha=0.05, i=2):
+        ci = []
+        for xi in np.atleast_1d(x).ravel():
+            try:
+                Lx = self.profile_quantile(xi, i=2)
+                ci.append(Lx.get_bounds(alpha=alpha))
+            except Exception:
+                ci.append((np.nan, np.nan))
+        return np.array(ci)
+
     def plotfitsummary(self, axes=None, fig=None):
         ''' Plot various diagnostic plots to asses the quality of the fit.
 
@@ -1403,20 +1423,10 @@ class FitDistribution(rv_frozen):
                       self.data, self.sf(self.data), symb1)
 
         if True:
-            i0 = 2
             low = int(np.log10(1.0/n)-0.7) - 1
-            log_sf = np.log(np.logspace(low, -0.5, 7)[::-1])
+            sf1 = np.logspace(low, -0.5, 7)[::-1]
-            T = self.isf(np.exp(log_sf))
+            ci1 = self.ci_sf(sf, alpha=0.05, i=2)
-            ci = []
+            axis.semilogy(self.isf(sf1), ci1, 'r--')
-            t = []
-            for Ti, log_sfi in zip(T, log_sf):
-                try:
-                    Lx = self.profile_probability(log_sfi, i=i0)
-                    ci.append(np.exp(Lx.get_bounds(alpha=0.05)))
-                    t.append(Ti)
-                except:
-                    pass
-            axis.semilogy(t, ci, 'r--')
         axis.set_xlabel('x')
         axis.set_ylabel('F(x) (%s)' % self.dist.name)
         axis.set_title('Empirical SF plot')

wafo/tests/test_misc.py

@@ -73,7 +73,6 @@ def test_disufq():
         # deep water
         assert_array_almost_equal(data, data2)
         assert_array_almost_equal(data, d_truth[i])
-
     # assert(False)