diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py
index baee24d..de82919 100644
--- a/wafo/spectrum/core.py
+++ b/wafo/spectrum/core.py
@@ -1422,8 +1422,14 @@ class SpecData1D(PlotData):
         The joint density of zero separated Max2min cycles in time (a);
         in space (b); AcAt in time for nonlinear sea model (c):
 
-        Hm0=7;Tp=11
-        S = jonswap(4*pi/Tp,[Hm0 Tp])
+        >>> from wafo.spectrum import models as sm
+        >>> w = np.linspace(0,3,100)
+        >>> Sj = sm.Jonswap()
+        >>> S = Sj.tospecdata()
+        >>> f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7)
+        >>> Hm0 = 7
+        >>> Tp = 11
+        >>> S = sm.Jonswap(4*pi/Tp,[Hm0 Tp])
         Sk = spec2spec(S,'k1d')
         L0 = spec2mom(S,1)
         paramu = [sqrt(L0)*[-4 4] 41]
diff --git a/wafo/stats/__init__.py b/wafo/stats/__init__.py
index e12a20a..0e5f112 100644
--- a/wafo/stats/__init__.py
+++ b/wafo/stats/__init__.py
@@ -341,6 +341,7 @@ from .core import *
 from .distributions import *
 from . import estimation
 
+
 # remove vonmises_cython from __all__, I don't know why it is included
 __all__ = [s for s in dir() if not (s.startswith('_') or s.endswith('cython'))]
 
diff --git a/wafo/stats/_distn_infrastructure.py b/wafo/stats/_distn_infrastructure.py
index 8bc802c..bbaa418 100644
--- a/wafo/stats/_distn_infrastructure.py
+++ b/wafo/stats/_distn_infrastructure.py
@@ -289,7 +289,7 @@ def nlogps(self, theta, x):
     if not self._argcheck(*args) or scale <= 0:
         return inf
     x = asarray((x - loc) / scale)
-    cond0 = (x <= self.a) | (self.b <= x)
+    cond0 = (x < self.a) | (self.b < x)
     Nbad = np.sum(cond0)
     if Nbad > 0:
         x = argsreduce(~cond0, x)[0]
@@ -328,11 +328,40 @@ def nlogps(self, theta, x):
     return T
 
 
+def _penalized_nnlf(self, theta, x):
+        ''' Return negative loglikelihood function,
+        i.e., - sum (log pdf(x, theta), axis=0)
+           where theta are the parameters (including loc and scale)
+        '''
+        try:
+            loc = theta[-2]
+            scale = theta[-1]
+            args = tuple(theta[:-2])
+        except IndexError:
+            raise ValueError("Not enough input arguments.")
+        if not self._argcheck(*args) or scale <= 0:
+            return inf
+        x = asarray((x-loc) / scale)
+
+        loginf = log(_XMAX)
+
+        if np.isneginf(self.a).all() and np.isinf(self.b).all():
+            Nbad = 0
+        else:
+            cond0 = (x < self.a) | (self.b < x)
+            Nbad = sum(cond0)
+            if Nbad > 0:
+                x = argsreduce(~cond0, x)[0]
+
+        N = len(x)
+        return self._nnlf(x, *args) + N*log(scale) + Nbad * 100.0 * loginf
+
+
 def _reduce_func(self, args, options):
     # First of all, convert fshapes params to fnum: eg for stats.beta,
     # shapes='a, b'. To fix `a`, can specify either `f1` or `fa`.
     # Convert the latter into the former.
-    kwds = options.copy()
+    kwds = options  # .copy()
     if self.shapes:
         shapes = self.shapes.replace(',', ' ').split()
         for j, s in enumerate(shapes):
@@ -574,6 +603,7 @@ rv_continuous.freeze = freeze
 rv_continuous.link = link
 rv_continuous._link = _link
 rv_continuous.nlogps = nlogps
+rv_continuous._penalized_nnlf = _penalized_nnlf
 rv_continuous._reduce_func = _reduce_func
 rv_continuous.fit = fit
 rv_continuous.fit2 = fit2
diff --git a/wafo/stats/estimation.py b/wafo/stats/estimation.py
index 6fbd769..2e92c13 100644
--- a/wafo/stats/estimation.py
+++ b/wafo/stats/estimation.py
@@ -134,7 +134,9 @@ class Profile(object):
     >>> sf_ci = Lsf.get_bounds(alpha=0.2)
     '''
 
-    def __init__(self, fit_dist, **kwds):
+    def __init__(self, fit_dist, method=None, **kwds):
+        if method is None:
+            method = fit_dist.method.lower()
 
         try:
             i0 = (1 - np.isfinite(fit_dist.par_fix)).argmax()
@@ -153,10 +155,10 @@ class Profile(object):
             [i0, 100, 0.05, None, None, None, None, None])]
 
         self.title = '%g%s CI' % (100 * (1.0 - self.alpha), '%')
-        if fit_dist.method.startswith('ml'):
+        if method.startswith('ml'):
             self.ylabel = self.ylabel + 'likelihood'
             Lmax = fit_dist.LLmax
-        elif fit_dist.method.startswith('mps'):
+        elif method.startswith('mps'):
             self.ylabel = self.ylabel + ' product spacing'
             Lmax = fit_dist.LPSmax
         else:
@@ -604,13 +606,13 @@ class FitDistribution(rv_frozen):
         self.dist = dist
         numargs = dist.numargs
 
-        self.method = method
+        self.method = method.lower()
         self.alpha = alpha
         self.par_fix = par_fix
         self.search = search
         self.copydata = copydata
 
-        if self.method.lower()[:].startswith('mps'):
+        if self.method.startswith('mps'):
             self._fitfun = self._nlogps
         else:
             self._fitfun = self._nnlf
@@ -797,7 +799,7 @@ class FitDistribution(rv_frozen):
             return np.inf
         dist = self.dist
         x = asarray((x - loc) / scale)
-        cond0 = (x <= dist.a) | (dist.b <= x)
+        cond0 = (x < dist.a) | (dist.b < x)
         Nbad = np.sum(cond0)
         if Nbad > 0:
             x = argsreduce(~cond0, x)[0]
@@ -816,8 +818,8 @@ class FitDistribution(rv_frozen):
         if any(tie):
             # TODO : implement this method for treating ties in data:
             # Assume measuring error is delta. Then compute
-            # yL = F(xi-delta,theta)
-            # yU = F(xi+delta,theta)
+            # yL = F(xi - delta, theta)
+            # yU = F(xi + delta, theta)
             # and replace
             # logDj = log((yU-yL)/(r-1)) for j = i+1,i+2,...i+r-1
 
@@ -1020,7 +1022,7 @@ class FitDistribution(rv_frozen):
                 fixstr = 'Fixed: phat[%s] = %s ' % (phatistr, phatvstr)
 
         infostr = 'Fit method: %s, Fit p-value: %2.2f %s' % (
-            self.method, self.pvalue, fixstr)
+            self.method.upper(), self.pvalue, fixstr)
         try:
             plotbackend.figtext(0.05, 0.01, infostr)
         except: