diff --git a/wafo/tests/test_gaussian.py b/wafo/tests/test_gaussian.py
index 3f7da42..d3faef8 100644
--- a/wafo/tests/test_gaussian.py
+++ b/wafo/tests/test_gaussian.py
@@ -5,70 +5,77 @@ Created on 17. juli 2010
 '''
 import numpy as np  # @UnusedImport
 from numpy import pi, inf  # @UnusedImport
-# @UnusedImport
+from numpy.testing import assert_array_almost_equal
 from wafo.gaussian import (Rind, prbnormtndpc, prbnormndpc, prbnormnd,
                            cdfnorm2d, prbnorm2d)
 
 
 def test_rind():
-    '''
-    >>> Et = 0.001946 # #  exact prob.
-    >>> n = 5
-    >>> Blo =-np.inf; Bup=-1.2; indI=[-1, n-1]  # Barriers
-    >>> m = np.zeros(n); rho = 0.3;
-    >>> Sc =(np.ones((n,n))-np.eye(n))*rho+np.eye(n)
-    >>> rind = Rind()
-    >>> E0, err0, terr0 = rind(Sc,m,Blo,Bup,indI);
-
-    >>> np.abs(E0-Et)< err0+terr0
-    array([ True], dtype=bool)
-    >>> 'E0 = %2.6f' % E0
-    'E0 = 0.001946'
-
-    >>> A = np.repeat(Blo,n); B = np.repeat(Bup,n)  # Integration limits
-    >>> E1, err1, terr1  = rind(np.triu(Sc),m,A,B); #same as E0
-    >>> np.abs(E1-Et)< err0+terr0
-    array([ True], dtype=bool)
-    >>> 'E1 = %2.5f' % E1
-    'E1 = 0.00195'
-
-    Compute expectation E( abs(X1*X2*...*X5) )
-    >>> xc = np.zeros((0,1))
-    >>> infinity = 37
-    >>> dev = np.sqrt(np.diag(Sc))  # std
-    >>> ind = np.nonzero(indI[1:])[0]
-    >>> Bup, Blo = np.atleast_2d(Bup,Blo)
-    >>> Bup[0,ind] = np.minimum(Bup[0,ind] , infinity*dev[indI[ind+1]])
-    >>> Blo[0,ind] = np.maximum(Blo[0,ind] ,-infinity*dev[indI[ind+1]])
-    >>> rind(Sc,m,Blo,Bup,indI, xc, nt=0)
-    (array([ 0.05494076]), array([ 0.00083066]), array([  1.00000000e-10]))
-
-    Compute expectation E( X1^{+}*X2^{+} ) with random
-    correlation coefficient,Cov(X1,X2) = rho2.
-    >>> m2  = [0, 0]
-    >>> rho2 = 0.3 #np.random.rand(1)
-    >>> Sc2 = [[1, rho2], [rho2 ,1]]
-    >>> Blo2 = 0; Bup2 = np.inf; indI2 = [-1, 1]
-    >>> rind2 = Rind(method=1)
-    >>> g2 = lambda x : (x*(np.pi/2+np.arcsin(x))+np.sqrt(1-x**2))/(2*np.pi)
-    >>> E2 = g2(rho2); E2   # exact value
-    0.24137214191774381
-
-    >>> E3, err3, terr3 = rind(Sc2,m2,Blo2,Bup2,indI2,nt=0); E3;err3;terr3
-    array([ 0.24127499])
-    array([ 0.00013838])
-    array([  1.00000000e-10])
-
-    >>> E4, err4, terr4 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0); E4;err4;terr4
-    array([ 0.24127499])
-    array([ 0.00013838])
-    array([  1.00000000e-10])
-
-    >>> E5, err5, terr5 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0,abseps=1e-4); E5;err5;terr5
-    array([ 0.24127499])
-    array([ 0.00013838])
-    array([  1.00000000e-10])
-    '''
+
+    Et = 0.001946  # exact prob.
+    n = 5
+    Blo = -np.inf
+    Bup = -1.2
+    indI = [-1, n - 1]  # Barriers
+    m = np.zeros(n)
+    rho = 0.3
+    Sc = (np.ones((n, n)) - np.eye(n)) * rho + np.eye(n)
+    rind = Rind()
+    E0, err0, terr0 = rind(Sc, m, Blo, Bup, indI)
+
+    assert(np.abs(E0 - Et) < err0 + terr0)
+
+    t = 'E0 = %2.6f' % E0
+    assert(t == 'E0 = 0.001946')
+
+    A = np.repeat(Blo, n)
+    B = np.repeat(Bup, n)  # Integration limits
+    E1, err1, terr1 = rind(np.triu(Sc), m, A, B)  # same as E0
+    assert(np.abs(E1 - Et) < err0 + terr0)
+
+    t = 'E1 = %2.5f' % E1
+    assert(t == 'E1 = 0.00195')
+
+    # Compute expectation E( abs(X1*X2*...*X5) )
+    xc = np.zeros((0, 1))
+    infinity = 37
+    dev = np.sqrt(np.diag(Sc))  # std
+    ind = np.nonzero(indI[1:])[0]
+    Bup, Blo = np.atleast_2d(Bup, Blo)
+    Bup[0, ind] = np.minimum(Bup[0, ind], infinity * dev[indI[ind + 1]])
+    Blo[0, ind] = np.maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]])
+    val, err, terr = rind(Sc, m, Blo, Bup, indI, xc, nt=0)
+    assert_array_almost_equal(val, 0.05494076)
+    assert(err < 0.001)
+    assert_array_almost_equal(terr, 1.00000000e-10)
+
+    # Compute expectation E( X1^{+}*X2^{+} ) with random
+    # correlation coefficient,Cov(X1,X2) = rho2.
+    m2 = [0, 0]
+    rho2 = 0.3  # np.random.rand(1)
+    Sc2 = [[1, rho2], [rho2, 1]]
+    Blo2 = 0
+    Bup2 = np.inf
+    indI2 = [-1, 1]
+    rind2 = Rind(method=1)
+    g2 = lambda x: (
+        x * (np.pi / 2 + np.arcsin(x)) + np.sqrt(1 - x**2)) / (2 * np.pi)
+    assert_array_almost_equal(g2(rho2), 0.24137214191774381)  # exact value
+
+    E3, err3, terr3 = rind(Sc2, m2, Blo2, Bup2, indI2, nt=0)
+    assert_array_almost_equal(E3, 0.24127499)
+    assert_array_almost_equal(err3, 0.00013838)
+    assert_array_almost_equal(terr3, 1.00000000e-10)
+
+    E4, err4, terr4 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0)
+    assert_array_almost_equal(E4, 0.24127499)
+    assert_array_almost_equal(err4, 0.00013838)
+    assert_array_almost_equal(terr4, 1.00000000e-10)
+#
+#     >>> E5, err5, terr5 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0,abseps=1e-4)
+#     array([ 0.24127499])
+#     array([ 0.00013838])
+#     array([  1.00000000e-10])
 
 
 def test_prbnormtndpc():
@@ -85,9 +92,9 @@ def test_prbnormtndpc():
     >>> rho3 = np.random.rand(3)
     >>> a3   = np.zeros(3)
     >>> b3   = np.repeat(inf,3)
-    >>> [val3,err3, ift3] = prbnormtndpc(rho3,a3,b3)  
+    >>> [val3,err3, ift3] = prbnormtndpc(rho3,a3,b3)
     >>> g3 = lambda x : 0.5-sum(np.sort(np.arccos([x[0]*x[1],x[0]*x[2],x[1]*x[2]])))/(4*pi)
-    >>> E3 = g3(rho3)   #  Exact value  
+    >>> E3 = g3(rho3)   #  Exact value
     >>> np.abs(E3-val3)<err3
     True
     '''