pep8 + some refactoring of XXXX.link() methods...

pywafo/src/wafo/stats/_distn_infrastructure.py

@@ -23,7 +23,7 @@ from scipy import optimize
 from scipy import integrate
 
 # to approximate the pdf of a continuous distribution given its cdf
-from scipy.misc import comb, derivative
+from scipy.misc import comb, derivative  # @UnresolvedImport
 
 from numpy import (arange, putmask, ravel, take, ones, sum, shape,
                    product, reshape, zeros, floor, logical_and, log, sqrt, exp,
@@ -1553,7 +1553,7 @@ class rv_continuous(rv_generic):
             for item in ['callparams', 'default', 'before_notes']:
                 tempdict[item] = tempdict[item].replace(
                     "\n%(shapes)s : array_like\n    shape parameters", "")
-        for i in range(2):
+        for _i in range(2):
             if self.shapes is None:
                 # necessary because we use %(shapes)s in two forms (w w/o ", ")
                 self.__doc__ = self.__doc__.replace("%(shapes)s, ", "")
@@ -1954,12 +1954,39 @@ class rv_continuous(rv_generic):
         return output
 
     def link(self, x, logSF, theta, i):
-        ''' Return dist. par. no. i as function of quantile (x) and log survival probability (sf)
-            where
-            theta is the list containing all parameters including location and scale.
         ''' 
-        raise ValueError('Link function not implemented for the %s distribution' % self.name)
-        return None
+        Return theta[i] as function of quantile, survival probability and
+            theta[j] for j!=i.
+
+        Parameters
+        ----------
+        x : quantile
+        logSF : logarithm of the survival probability
+        theta : list
+            all distribution parameters including location and scale.
+
+        Returns
+        -------
+        theta[i] : real scalar
+            fixed distribution parameter theta[i] as function of x, logSF and
+            theta[j] where j != i.
+
+        LINK is a function connecting the fixed distribution parameter theta[i]
+        with the quantile (x) and the survival probability (SF) and the
+        remaining free distribution parameters theta[j] for j!=i, i.e.:
+            theta[i] = link(x, logSF, theta, i),
+        where logSF = log(Prob(X>x; theta)).
+
+        See also 
+        estimation.Profile
+        '''
+        return self._link(x, logSF, theta, i)
+
+    def _link(self, x, logSF, theta, i):
+        msg = ('Link function not implemented for the %s distribution' %
+               self.name)
+        raise NotImplementedError(msg)
+
 
     def _nnlf(self, x, *args):
         return -sum(self._logpdf(x, *args), axis=0)
@@ -3030,7 +3057,7 @@ class rv_discrete(rv_generic):
             for item in ['callparams', 'default', 'before_notes']:
                 tempdict[item] = tempdict[item].replace(
                     "\n%(shapes)s : array_like\n    shape parameters", "")
-        for i in range(2):
+        for _i in range(2):
             if self.shapes is None:
                 # necessary because we use %(shapes)s in two forms (w w/o ", ")
                 self.__doc__ = self.__doc__.replace("%(shapes)s, ", "")
@@ -3498,7 +3525,7 @@ class rv_discrete(rv_generic):
         else:
             invfac = 1.0
 
-        tot = 0.0
+        #tot = 0.0
         low, upp = self._ppf(0.001, *args), self._ppf(0.999, *args)
         low = max(min(-suppnmin, low), lb)
         upp = min(max(suppnmin, upp), ub)

pywafo/src/wafo/stats/estimation.py

@@ -152,8 +152,8 @@ class Profile(object):
         with ML or MPS estimated distribution parameters.
     **kwds : named arguments with keys
         i : scalar integer
-            defining which distribution parameter to profile, i.e. which
-            parameter to keep fixed (default first non-fixed parameter)
+            defining which distribution parameter to keep fixed in the 
+            profiling process (default first non-fixed parameter)
         pmin, pmax : real scalars
             Interval for either the parameter, phat(i), prb, or x, used in the
             optimization of the profile function (default is based on the
@@ -236,22 +236,21 @@ class Profile(object):
         self.fit_dist = fit_dist
         self.data = None
         self.args = None
-        self.title = 'Profile log'
+        self.title = ''
         self.xlabel = ''
-        self.ylabel = ''
+        self.ylabel = 'Profile log'
         (self.i_fixed, self.N, self.alpha, self.pmin, self.pmax, self.x,
          self.logSF, self.link) = map(
             kwds.get,
-            ['i', 'N', 'alpha', 'pmin',
-                                'pmax', 'x', 'logSF', 'link'],
+            ['i', 'N', 'alpha', 'pmin', 'pmax', 'x', 'logSF', 'link'],
             [i0, 100, 0.05, None, None, None, None, None])
 
-        self.ylabel = '%g%s CI' % (100 * (1.0 - self.alpha), '%')
+        self.title = '%g%s CI' % (100 * (1.0 - self.alpha), '%')
         if fit_dist.method.startswith('ml'):
-            self.title = self.title + 'likelihood'
+            self.ylabel = self.ylabel + 'likelihood'
             Lmax = fit_dist.LLmax
         elif fit_dist.method.startswith('mps'):
-            self.title = self.title + ' product spacing'
+            self.ylabel = self.ylabel + ' product spacing'
             Lmax = fit_dist.LPSmax
         else:
             raise ValueError(
@@ -267,7 +266,7 @@ class Profile(object):
         self.i_fixed = atleast_1d(self.i_fixed)
 
         if 1 - isnotfixed[self.i_fixed]:
-            raise ValueError(
+            raise IndexError(
                 "Index i must be equal to an index to one of the free " +
                 "parameters.")
 
@@ -295,7 +294,7 @@ class Profile(object):
             self.xlabel = 'phat(%d)' % self.i_fixed
             p_opt = self._par[self.i_fixed]
         elif self.profile_x:
-            self.logSF = log(fit_dist.sf(self.x))
+            self.logSF = fit_dist.logsf(self.x)
             self._local_link = lambda fix_par, par: self.link(
                 fix_par, self.logSF, par, self.i_fixed)
             self.xlabel = 'x'
@@ -522,16 +521,23 @@ class Profile(object):
             CI = ecross(self.args, self.data, ind[[0, -1]], cross_level)
         return CI
 
-    def plot(self):
+    def plot(self, axis=None):
         ''' Plot profile function with 100(1-alpha)% CI
         '''
-        plotbackend.plot(
+        if axis is None:
+            axis = plotbackend.gca()
+            
+        p_ci = self.get_bounds(self.alpha)
+        axis.plot(
             self.args, self.data,
-            self.args[[0, -1]], [self.Lmax, ] * 2, 'r',
-            self.args[[0, -1]], [self.alpha_cross_level, ] * 2, 'r')
-        plotbackend.title(self.title)
-        plotbackend.ylabel(self.ylabel)
-        plotbackend.xlabel(self.xlabel)
+            self.args[[0, -1]], [self.Lmax, ] * 2, 'r--',
+            self.args[[0, -1]], [self.alpha_cross_level, ] * 2, 'r--')
+        axis.vlines(p_ci, ymin=axis.get_ylim()[0],
+                           ymax=self.Lmax, #self.alpha_cross_level,
+                           color='r', linestyles='--')
+        axis.set_title(self.title)
+        axis.set_ylabel(self.ylabel)
+        axis.set_xlabel(self.xlabel)
 
 
 def _discretize_adaptive(fun, a, b, tol=0.005, n=5):