diff --git a/pywafo/src/wafo/stats/distributions.py b/pywafo/src/wafo/stats/distributions.py
index 02040ef..49db88c 100644
--- a/pywafo/src/wafo/stats/distributions.py
+++ b/pywafo/src/wafo/stats/distributions.py
@@ -2031,6 +2031,64 @@ class rv_continuous(rv_generic):
         else:
             N = len(x)
             return self._nnlf(x, *args) + N*log(scale)
+    def hessian_nlogps(self, theta, data, eps=None):
+        ''' approximate hessian of nlogps where theta are the parameters (including loc and scale)
+        '''
+        #Nd = len(x)
+        np = len(theta)
+        # pab 07.01.2001: Always choose the stepsize h so that
+        # it is an exactly representable number.
+        # This is important when calculating numerical derivatives and is
+        #  accomplished by the following.
+
+        if eps == None:
+            eps = (floatinfo.machar.eps) ** 0.4
+        #xmin = floatinfo.machar.xmin
+        #myfun = lambda y: max(y,100.0*log(xmin)) #% trick to avoid log of zero
+        delta = (eps + 2.0) - 2.0
+        delta2 = delta ** 2.0
+        #    % Approximate 1/(nE( (d L(x|theta)/dtheta)^2)) with
+        #    %             1/(d^2 L(theta|x)/dtheta^2)
+        #    %  using central differences
+
+        LL = self.nlogps(theta, data)
+        H = zeros((np, np))   #%% Hessian matrix
+        theta = tuple(theta)
+        for ix in xrange(np):
+            sparam = list(theta)
+            sparam[ix] = theta[ix] + delta
+            fp = self.nlogps(sparam, data)
+            #fp = sum(myfun(x))
+
+            sparam[ix] = theta[ix] - delta
+            fm = self.nlogps(sparam, data)
+            #fm = sum(myfun(x))
+
+            H[ix, ix] = (fp - 2 * LL + fm) / delta2
+            for iy in range(ix + 1, np):
+                sparam[ix] = theta[ix] + delta
+                sparam[iy] = theta[iy] + delta
+                fpp = self.nlogps(sparam, data)
+                #fpp = sum(myfun(x))
+
+                sparam[iy] = theta[iy] - delta
+                fpm = self.nlogps(sparam, data)
+                #fpm = sum(myfun(x))
+
+                sparam[ix] = theta[ix] - delta
+                fmm = self.nlogps(sparam, data)
+                #fmm = sum(myfun(x));
+
+                sparam[iy] = theta[iy] + delta
+                fmp = self.nlogps(sparam, data)
+                #fmp = sum(myfun(x))
+                H[ix, iy] = ((fpp + fmm) - (fmp + fpm)) / (4. * delta2)
+                H[iy, ix] = H[ix, iy]
+                sparam[iy] = theta[iy];
+
+        # invert the Hessian matrix (i.e. invert the observed information number)
+        #pcov = -pinv(H);
+        return -H    
     def hessian_nnlf(self, theta, data, eps=None):
         ''' approximate hessian of nnlf where theta are the parameters (including loc and scale)
         '''
diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py
index 075b24f..0983a66 100644
--- a/pywafo/src/wafo/stats/estimation.py
+++ b/pywafo/src/wafo/stats/estimation.py
@@ -192,7 +192,6 @@ class Profile(object):
     Examples
     --------
     # MLE 
-    >>> import numpy as np
     >>> import wafo.stats as ws
     >>> R = ws.weibull_min.rvs(1,size=100);
     >>> phat = FitDistribution(ws.weibull_min, R, 1, scale=1, floc=0.0)
@@ -364,6 +363,10 @@ class Profile(object):
                 pcov = self.fit_dist.par_cov[i_notfixed, :][:, i_notfixed]
                 pvar = sum(numpy.dot(drl, pcov) * drl)
 
+            if pvar<0: 
+                pvar = -pvar*2
+            elif numpy.isnan(pvar):
+                pvar = p_opt*0.5
             p_crit = -norm_ppf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5
             if self.pmin == None:
                 self.pmin = p_opt - 5.0 * p_crit
@@ -515,7 +518,6 @@ class FitDistribution(rv_frozen):
     --------
     Estimate distribution parameters for weibull_min distribution.
     >>> import wafo.stats as ws
-    >>> import numpy as np
     >>> R = ws.weibull_min.rvs(1,size=100);
     >>> phat = FitDistribution(ws.weibull_min, R, 1, scale=1, floc=0.0)
     
@@ -719,7 +721,8 @@ class FitDistribution(rv_frozen):
         '''Compute covariance
         '''
         somefixed = (self.par_fix != None) and any(isfinite(self.par_fix))
-        H = numpy.asmatrix(self.dist.hessian_nnlf(self.par, self.data))
+        H1 = numpy.asmatrix(self.dist.hessian_nnlf(self.par, self.data))
+        H = numpy.asmatrix(self.dist.hessian_nlogps(self.par, self.data))
         self.H = H
         try:
             if somefixed:
@@ -798,7 +801,6 @@ class FitDistribution(rv_frozen):
         Examples
         --------
         # MLE 
-        >>> import numpy as np
         >>> import wafo.stats as ws
         >>> R = ws.weibull_min.rvs(1,size=100);
         >>> phat = FitDistribution(ws.weibull_min, R, 1, scale=1, floc=0.0)
@@ -1006,13 +1008,15 @@ class FitDistribution(rv_frozen):
  
 
 def main():
-     
-#    import wafo.stats as ws
-#    R = ws.weibull_min.rvs(1,size=100);
-#    phat = FitDistribution(ws.weibull_min, R, 1, scale=1, floc=0.0)
-#    
+    import doctest
+    doctest.testmod()
+def test1():
+    import wafo.stats as ws
+    R = ws.weibull_min.rvs(1,size=100);
+    phat = FitDistribution(ws.weibull_min, R, method='mps')
+    
 #    # Better CI for phat.par[i=0]
-#    Lp1 = Profile(phat, i=0)
+    Lp1 = Profile(phat, i=1)
 #    Lp2 = Profile(phat, i=2)
 #    SF = 1./990
 #    x = phat.isf(SF)
@@ -1029,8 +1033,7 @@ def main():
 #    pass
     
     
-    import doctest
-    doctest.testmod()
+ 
     
     
 #    _WAFODIST = ppimport('wafo.stats.distributions')
@@ -1080,4 +1083,5 @@ def main():
 #    lp = pht.profile()
                 
 if __name__ == '__main__':
+    test1()
     main()