From deef6c4994bf088fb4a44b533ea1f27cbc4598a4 Mon Sep 17 00:00:00 2001
From: Per A Brodtkorb <per.andreas.brodtkorb@gmail.com>
Date: Fri, 14 Oct 2016 08:48:38 +0200
Subject: [PATCH] Added doctest to check_random_state and added ci_quantile and
 ci_sf to FitDistribution

---
 wafo/misc.py             | 22 ++++++++++++++++------
 wafo/stats/estimation.py | 36 +++++++++++++++++++++++-------------
 wafo/tests/test_misc.py  |  1 -
 3 files changed, 39 insertions(+), 20 deletions(-)

diff --git a/wafo/misc.py b/wafo/misc.py
index 859c116..2ea77b0 100644
--- a/wafo/misc.py
+++ b/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'
-                     ' instance' % seed)
+    msg = '{} cannot be used to seed a numpy.random.RandomState instance'
+    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,
-                                       linalg.lstsq(I-np.dot(Bp, Ap), e)[0])
+                        rh = I - np.dot(Bp, Ap)
+                        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__":
diff --git a/wafo/stats/estimation.py b/wafo/stats/estimation.py
index bfc8a3b..b94df59 100644
--- a/wafo/stats/estimation.py
+++ b/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])
-            T = self.isf(np.exp(log_sf))
-            ci = []
-            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--')
+            sf1 = np.logspace(low, -0.5, 7)[::-1]
+            ci1 = self.ci_sf(sf, alpha=0.05, i=2)
+            axis.semilogy(self.isf(sf1), ci1, 'r--')
         axis.set_xlabel('x')
         axis.set_ylabel('F(x) (%s)' % self.dist.name)
         axis.set_title('Empirical SF plot')
diff --git a/wafo/tests/test_misc.py b/wafo/tests/test_misc.py
index a234b7d..b6b1424 100644
--- a/wafo/tests/test_misc.py
+++ b/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)