diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
index 7559fe7..62a4e63 100644
--- a/.github/workflows/run_tests.yml
+++ b/.github/workflows/run_tests.yml
@@ -7,6 +7,11 @@ on:
   pull_request:
     branches:
       - main
+      - development
+    paths:
+      - "tutorials/**"
+      - "dlordinal/**"
+      - ".github/workflows/**"
 
 jobs:
   tests:
@@ -14,10 +19,10 @@ jobs:
 
     steps:
     - name: Checkout repository
-      uses: actions/checkout@v2
+      uses: actions/checkout@v4
 
     - name: Set up Python
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v5
       with:
         python-version: 3.8
 
@@ -38,10 +43,10 @@ jobs:
 
     steps:
     - name: Checkout repository
-      uses: actions/checkout@v2
+      uses: actions/checkout@v4
 
     - name: Set up Python
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v5
       with:
         python-version: 3.8
 
@@ -55,6 +60,7 @@ jobs:
     - name: Run tests for codecov
       run: |
         pytest --cov=dlordinal --cov-report=xml
+      timeout-minutes: 20
 
     - name: Upload coverage reports to Codecov
       uses: codecov/codecov-action@v3
diff --git a/dlordinal/datasets/tests/test_adience.py b/dlordinal/datasets/tests/test_adience.py
index ab281d1..a1388a6 100644
--- a/dlordinal/datasets/tests/test_adience.py
+++ b/dlordinal/datasets/tests/test_adience.py
@@ -8,7 +8,7 @@
 import pytest
 from PIL import Image
 
-from ..adience import Adience
+from dlordinal.datasets import Adience
 
 temp_dir = None
 
diff --git a/dlordinal/datasets/tests/test_featuredataset.py b/dlordinal/datasets/tests/test_featuredataset.py
index 591b540..2642688 100644
--- a/dlordinal/datasets/tests/test_featuredataset.py
+++ b/dlordinal/datasets/tests/test_featuredataset.py
@@ -2,7 +2,7 @@
 import pytest
 import torch
 
-from ..featuredataset import FeatureDataset
+from dlordinal.datasets import FeatureDataset
 
 
 @pytest.fixture
diff --git a/dlordinal/datasets/tests/test_fgnet.py b/dlordinal/datasets/tests/test_fgnet.py
index 3b1281c..cd073cb 100644
--- a/dlordinal/datasets/tests/test_fgnet.py
+++ b/dlordinal/datasets/tests/test_fgnet.py
@@ -3,7 +3,7 @@
 
 import pytest
 
-from ..fgnet import FGNet
+from dlordinal.datasets import FGNet
 
 TMP_DIR = "./tmp_test_dir_fgnet"
 
diff --git a/dlordinal/distributions/__init__.py b/dlordinal/distributions/__init__.py
index 439ad21..571772a 100644
--- a/dlordinal/distributions/__init__.py
+++ b/dlordinal/distributions/__init__.py
@@ -1,15 +1,19 @@
-from .beta_distribution import get_beta_probabilities
-from .binomial_distribution import get_binomial_probabilities
-from .exponential_distribution import get_exponential_probabilities
-from .general_triangular_distribution import get_general_triangular_probabilities
+from .beta_distribution import get_beta_softlabels
+from .binomial_distribution import get_binomial_softlabels
+from .exponential_distribution import get_exponential_softlabels
+from .general_triangular_distribution import (
+    get_general_triangular_params,
+    get_general_triangular_softlabels,
+)
 from .poisson_distribution import get_poisson_probabilities
-from .triangular_distribution import get_triangular_probabilities
+from .triangular_distribution import get_triangular_softlabels
 
 __all__ = [
-    "get_beta_probabilities",
-    "get_exponential_probabilities",
-    "get_binomial_probabilities",
+    "get_beta_softlabels",
+    "get_exponential_softlabels",
+    "get_binomial_softlabels",
     "get_poisson_probabilities",
-    "get_triangular_probabilities",
-    "get_general_triangular_probabilities",
+    "get_triangular_softlabels",
+    "get_general_triangular_params",
+    "get_general_triangular_softlabels",
 ]
diff --git a/dlordinal/distributions/beta_distribution.py b/dlordinal/distributions/beta_distribution.py
index a2f736a..3bb7694 100644
--- a/dlordinal/distributions/beta_distribution.py
+++ b/dlordinal/distributions/beta_distribution.py
@@ -40,7 +40,7 @@ def beta_dist(x, p, q, a=1.0):
     q : float
             Second shape parameter (:math:`q > 0`).
     a : float, default=1.0
-            Scaling parameter.
+            Scaling parameter (:math:`a > 0`).
 
     Returns
     -------
@@ -57,27 +57,60 @@ def beta_dist(x, p, q, a=1.0):
     return (x ** (a * p)) / (p * beta_func(p, q)) * hyp2f1(p, 1 - q, p + 1, x**a)
 
 
-def get_beta_probabilities(n, p, q, a=1.0):
-    """Get probabilities from a beta distribution :math:`B(p,q,a)` for ``n`` splits.
+def get_beta_softlabels(J, p, q, a=1.0):
+    """Get probabilities from a beta distribution :math:`B(p,q,a)` for ``J`` splits.
+    The :math:`[0,1]` interval is split into ``J`` intervals and the probability for
+    each interval is computed as the difference between the value of the distribution
+    function in the upper limit of the interval and the value of the distribution
+    function in the lower limit of the interval. Thus, the probability for class ``j``
+    is computed as :math:`B(p,q,a)(j/J) - B(p,q,a)((j-1)/J)`.
 
     Parameters
     ----------
-    n : int
-            Number of classes.
+    J : int
+            Number of classes or splits.
     p : float
             First shape parameter (:math:`p > 0`).
     q : float
             Second shape parameter (:math:`q > 0`).
     a : float, default=1.0
-            Scaling parameter.
+            Scaling parameter (:math:`a > 0`).
+
+    Raises
+    ------
+    ValueError
+            If ``J`` is not a positive integer, if ``p`` is not positive, if ``q`` is
+            not positive or if ``a`` is not positive.
 
     Returns
     -------
     probs: list
-            List of probabilities.
+            List of ``J`` elements that represent the probability associated with each
+            class or split.
+
+    Example
+    -------
+    >>> from dlordinal.distributions import get_beta_probabilities
+    >>> get_beta_probabilities(3, 2, 3)
+    [0.4074074080000002, 0.48148148059259255, 0.11111111140740726]
+    >>> get_beta_probabilities(5, 5, 1, a=2)
+    [1.0240000307200007e-07, 0.00010475520052121611, 0.005941759979320316,
+      0.10132756401484902, 0.8926258084053067]
     """
 
-    intervals = get_intervals(n)
+    if J < 2 or not isinstance(J, int):
+        raise ValueError(f"{J=} must be a positive integer greater than 1")
+
+    if p <= 0:
+        raise ValueError(f"{p=} must be positive")
+
+    if q <= 0:
+        raise ValueError(f"{q=} must be positive")
+
+    if a <= 0:
+        raise ValueError(f"{a=} must be positive")
+
+    intervals = get_intervals(J)
     probs = []
 
     # Compute probability for each interval (class) using the distribution function.
diff --git a/dlordinal/distributions/binomial_distribution.py b/dlordinal/distributions/binomial_distribution.py
index 0196712..893d4e1 100644
--- a/dlordinal/distributions/binomial_distribution.py
+++ b/dlordinal/distributions/binomial_distribution.py
@@ -2,22 +2,39 @@
 from scipy.stats import binom
 
 
-def get_binomial_probabilities(n):
-    """Get probabilities from binominal distribution for n classes.
+def get_binomial_softlabels(J):
+    """Get probabilities for the binomial distribution for ``J`` classes or splits
+    using the approach described in :footcite:t:`liu2020unimodal`.
+    The :math:`[0,1]` interval is split into ``J`` intervals and the probability for
+    each interval is computed as the difference between the value of the binomial
+    probability function for the interval boundaries. The probability for the first
+    interval is computed as the value of the binomial probability function for the first
+    interval boundary.
+
+    The binomial distributions employed are denoted as :math:`\\text{b}(k, n-1, p)` where
+    :math:`k` is given by the order of the class for which the probability is computed,
+    and :math:`p` is given by :math:`0.1 + (0.9-0.1) / (n-1) * j` where :math:`j` is
+    is the order of the target class.
 
     Parameters
     ----------
-    n : int
-            Number of classes.
+    J : int
+            Number of classes or splits.
+
+    Raises
+    ------
+    ValueError
+            If ``J`` is not a positive integer greater than 1.
 
     Returns
     -------
-    probs : 2d array-like
+    probs : 2d array-like of shape (J, J)
             Matrix of probabilities where each row represents the true class
-            and each column the probability for class n.
+            and each column the probability for class ``j``.
 
     Example
     -------
+    >>> from dlordinal.distributions import get_binominal_probabilities
     >>> get_binominal_probabilities(5)
     array([[6.561e-01, 2.916e-01, 4.860e-02, 3.600e-03, 1.000e-04],
             [2.401e-01, 4.116e-01, 2.646e-01, 7.560e-02, 8.100e-03],
@@ -26,6 +43,9 @@ def get_binomial_probabilities(n):
             [1.000e-04, 3.600e-03, 4.860e-02, 2.916e-01, 6.561e-01]])
     """
 
+    if J < 2 or not isinstance(J, int):
+        raise ValueError(f"{J=} must be a positive integer greater than 1")
+
     params = {}
 
     params["4"] = np.linspace(0.1, 0.9, 4)
@@ -39,7 +59,7 @@ def get_binomial_probabilities(n):
 
     probs = []
 
-    for true_class in range(0, n):
-        probs.append(binom.pmf(np.arange(0, n), n - 1, params[str(n)][true_class]))
+    for true_class in range(0, J):
+        probs.append(binom.pmf(np.arange(0, J), J - 1, params[str(J)][true_class]))
 
     return np.array(probs)
diff --git a/dlordinal/distributions/exponential_distribution.py b/dlordinal/distributions/exponential_distribution.py
index 738d2a7..c000431 100644
--- a/dlordinal/distributions/exponential_distribution.py
+++ b/dlordinal/distributions/exponential_distribution.py
@@ -2,26 +2,44 @@
 from scipy.special import softmax
 
 
-def get_exponential_probabilities(n, p=1.0, tau=1.0):
-    """Get probabilities from exponential distribution for n classes.
+def get_exponential_softlabels(J, p=1.0, tau=1.0):
+    """Get probabilities from exponential distribution for ``J`` classes or splits as
+    described in :footcite:t:`liu2020unimodal` and :footcite:t:`vargas2023exponential`.
+    The :math:`[0,1]` interval is split into ``J`` intervals and the probability for
+    each interval is computed as the difference between the value of the exponential
+    function for the interval boundaries. The probability for the first interval is
+    computed as the value of the exponential function for the first interval boundary.
+    Then, a softmax function is applied to each row of the resulting matrix to obtain
+    valid probabilities.
+
+    The aforementioned exponential function is defined as follows:
+    :math:`f(x; k, p, \\tau) = \\frac{-|x - k|^p}{\\tau}` where :math:`k` is given by the order of the
+    true class for which the probability is computed, :math:`p` is the exponent
+    parameter and :math:`\\tau` is the scaling parameter.
 
     Parameters
     ----------
-    n : int
+    J : int
             Number of classes.
     p : float, default=1.0
-            Exponent parameter.
+            Exponent parameter :math:`p`.
     tau : float, default=1.0
-            Scaling parameter.
+            Scaling parameter :math:`\\tau`.
+
+    Raises
+    ------
+    ValueError
+            If ``J`` is not a positive integer greater than 1.
 
     Returns
     -------
-    probs : 2d array-like
+    probs : 2d array-like of shape (J, J)
             Matrix of probabilities where each row represents the true class
-            and each column the probability for class n.
+            and each column the probability for class j.
 
     Example
     -------
+    >>> from dlordinal.distributions import get_exponential_probabilities
     >>> get_exponential_probabilities(5)
     array([[0.63640865, 0.23412166, 0.08612854, 0.03168492, 0.01165623],
     [0.19151597, 0.52059439, 0.19151597, 0.07045479, 0.02591887],
@@ -30,9 +48,12 @@ def get_exponential_probabilities(n, p=1.0, tau=1.0):
     [0.01165623, 0.03168492, 0.08612854, 0.23412166, 0.63640865]])
     """
 
+    if J < 2 or not isinstance(J, int):
+        raise ValueError(f"{J=} must be a positive integer greater than 1")
+
     probs = []
 
-    for true_class in range(0, n):
-        probs.append(-(np.abs(np.arange(0, n) - true_class) ** p) / tau)
+    for true_class in range(0, J):
+        probs.append(-(np.abs(np.arange(0, J) - true_class) ** p) / tau)
 
     return softmax(np.array(probs), axis=1)
diff --git a/dlordinal/distributions/general_triangular_distribution.py b/dlordinal/distributions/general_triangular_distribution.py
index 1ae1fbb..b5acfec 100644
--- a/dlordinal/distributions/general_triangular_distribution.py
+++ b/dlordinal/distributions/general_triangular_distribution.py
@@ -7,62 +7,98 @@
 from .utils import get_intervals, triangular_cdf
 
 
-def get_general_triangular_params(n: int, alphas: np.ndarray, verbose: int = 0):
+def get_general_triangular_params(J: int, alphas: np.ndarray, verbose: int = 0):
     """
-    Get the parameters for the general triangular distribution.
+    Get the parameters (:math:`a`, :math:`b` and :math:`c`) for the
+    general triangular distribution. The parameters are computed using
+    the :math:`\\boldsymbol{\\alpha}` values and the number of classes
+    :math:`J`. The :math:`\\boldsymbol{\\alpha}` vector contains two :math:`\\alpha`
+    values for each class or split. :math:`\\alpha_0` should always be 0, given that
+    the error on the left side of the first class is always 0. In the same way, the
+    :math:`\\alpha_{2J}` value should always be 0, given that the error on the right
+    side of the last class is always 0. The :math:`\\alpha` values for the other
+    classes should be between 0 and 1.
+
+    The parameters :math:`a`, :math:`b` and :math:`c` for class :math:`j` are computed
+    as described in :footcite:t:`vargas2023gentri`.
 
     Parameters
     ----------
-    n : int
-        Number of classes.
+    J : int
+        Number of classes or splits.
     alphas : np.ndarray
-        Array of alphas.
+        Array that represents the error on the left and right side of each class.
+        It is the :math:`\\boldsymbol{\\alpha}` vector described
+        in :footcite:t:`vargas2023gentri`.
     verbose : int, optional
         Verbosity level, by default 0.
+
+    Raises
+    ------
+    ValueError
+        If the number of classes :math:`J` is less than 2.
+        If the :math:`\\boldsymbol{\\alpha}` vector is not a numpy array of shape
+        :math:`(2J,)`.
+
+    Returns
+    -------
+    list
+        List of dictionaries with the parameters for each class.
+
+    Example
+    -------
+    >>> from dlordinal.distributions import get_general_triangular_params
+    >>> get_general_triangular_params(5, [0, 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.1, 0])
+    [{'alpha2j_1': 0, 'alpha2j': 0.05, 'a': 0, 'b': 0.25760143110525874, 'c': 0}, {'alpha2j_1': 0.05, 'alpha2j': 0.05, 'a': 0.153752470442574, 'b': 0.446247529557426, 'c': 0.3}, {'alpha2j_1': 0.05, 'alpha2j': 0.05, 'a': 0.353752470442574, 'b': 0.646247529557426, 'c': 0.5}, {'alpha2j_1': 0.05, 'alpha2j': 0.1, 'a': 0.550779686438060, 'b': 0.875486049708105, 'c': 0.7}, {'alpha2j_1': 0.0, 'alpha2j': 0, 'a': 0.8, 'b': 1, 'c': 1}]
+
     """
+
     alphas = np.array(alphas)
-    if not isinstance(alphas, np.ndarray) or alphas.shape != (2 * n,):
+    if not isinstance(alphas, np.ndarray) or alphas.shape != (2 * J,):
         raise ValueError(
             f"alphas must be a numpy array of shape (2 * n,), but got {alphas.shape}"
         )
 
-    def abc1(n, alpha):
+    if J < 2:
+        raise ValueError(f"J must be greater than 1, but got {J}")
+
+    def abc1(J, alpha):
         a = 0
-        b = 1.0 / ((1.0 - math.sqrt(alpha)) * n)
+        b = 1.0 / ((1.0 - math.sqrt(alpha)) * J)
         c = 0
         return a, b, c
 
-    def abcJ(n, alpha):
-        a = 1.0 - 1.0 / (n * (1.0 - math.sqrt(alpha)))
+    def abcJ(J, alpha):
+        a = 1.0 - 1.0 / (J * (1.0 - math.sqrt(alpha)))
         b = 1
         c = 1
         return a, b, c
 
-    def abcj(n, j, alpha2j_1, alpha2j):
+    def abcj(J, j, alpha2j_1, alpha2j):
         # Calculation of the terms of the system of equations
         # b = ...
-        numa2 = 2 * n**2 - 2 * n**2 * alpha2j_1
-        numa1 = 2 * j * n * alpha2j_1 - n * alpha2j_1 - 4 * n * j + 4 * n
+        numa2 = 2 * J**2 - 2 * J**2 * alpha2j_1
+        numa1 = 2 * j * J * alpha2j_1 - J * alpha2j_1 - 4 * J * j + 4 * J
         numa0 = 2 * j**2 - 4 * j + 2
-        dena0 = 2 * j * n * alpha2j_1 - n * alpha2j_1
-        dena1 = -(2 * n**2 * alpha2j_1)
+        dena0 = 2 * j * J * alpha2j_1 - J * alpha2j_1
+        dena1 = -(2 * J**2 * alpha2j_1)
 
         # a = ...
-        numb2 = 2 * n**2 * alpha2j - 2 * n**2
-        numb1 = 4 * n * j - 2 * n * j * alpha2j + n * alpha2j
+        numb2 = 2 * J**2 * alpha2j - 2 * J**2
+        numb1 = 4 * J * j - 2 * J * j * alpha2j + J * alpha2j
         numb0 = -(2 * j**2)
-        denb1 = 2 * n**2 * alpha2j
-        denb0 = -2 * n * j * alpha2j + n * alpha2j
+        denb1 = 2 * J**2 * alpha2j
+        denb0 = -2 * J * j * alpha2j + J * alpha2j
 
         a, b = symbols("a, b")
-        c = (2 * j - 1) / (2 * n)
+        c = (2 * j - 1) / (2 * J)
         eq1 = Eq((numa2 * a**2 + numa1 * a + numa0) / (dena0 + dena1 * a), b)
         eq2 = Eq((numb2 * b**2 + numb1 * b + numb0) / (denb1 * b + denb0), a)
 
         try:
             nsol = nsolve([eq1, eq2], [a, b], [1, 1])
 
-            if nsol[0] < (j - 1) / n and nsol[1] > j / n:
+            if nsol[0] < (j - 1) / J and nsol[1] > j / J:
                 return nsol[0], nsol[1], c
         except ValueError:
             pass
@@ -86,31 +122,31 @@ def abcj(n, j, alpha2j_1, alpha2j):
             if isinstance(s_b, Add):
                 s_b = s_b.args[0]
 
-            if s_a < (j - 1) / n and s_b > j / n:
+            if s_a < (j - 1) / J and s_b > j / J:
                 return s_a, s_b, c
 
         raise ValueError(f"Could not find solution for {j=}, {alpha2j_1=}, {alpha2j=}")
 
     if verbose >= 3:
-        print(f"{abc1(n,alphas[2])=}, {abcJ(n,alphas[2*n-1])=}")
-        for i in range(2, n):
-            print(f"{i=}  {abcj(n,i,alphas[2*i-1],alphas[2*i])}")
+        print(f"{abc1(J,alphas[2])=}, {abcJ(J,alphas[2*J-1])=}")
+        for i in range(2, J):
+            print(f"{i=}  {abcj(J,i,alphas[2*i-1],alphas[2*i])}")
 
     params = []
 
     # Compute params for each interval (class)
-    for j in range(1, n + 1):
+    for j in range(1, J + 1):
         j_params = {}
         if j == 1:
-            a, b, c = abc1(n, alphas[2])
+            a, b, c = abc1(J, alphas[2])
             j_params["alpha2j_1"] = 0
             j_params["alpha2j"] = alphas[2]
-        elif j == n:
-            a, b, c = abcJ(n, alphas[2 * n - 1])
-            j_params["alpha2j_1"] = alphas[2 * n - 1]
+        elif j == J:
+            a, b, c = abcJ(J, alphas[2 * J - 1])
+            j_params["alpha2j_1"] = alphas[2 * J - 1]
             j_params["alpha2j"] = 0
         else:
-            a, b, c = abcj(n, j, alphas[2 * j - 1], alphas[2 * j])
+            a, b, c = abcj(J, j, alphas[2 * j - 1], alphas[2 * j])
             j_params["alpha2j_1"] = alphas[2 * j - 1]
             j_params["alpha2j"] = alphas[2 * j]
 
@@ -122,25 +158,70 @@ def abcj(n, j, alpha2j_1, alpha2j):
     return params
 
 
-def get_general_triangular_probabilities(n: int, alphas: np.ndarray, verbose: int = 0):
-    """
-    Get the probabilities for the general triangular distribution.
+def get_general_triangular_softlabels(J: int, alphas: np.ndarray, verbose: int = 0):
+    """Get probabilities from triangular distributions for ``J`` classes or splits.
+    The :math:`[0,1]` interval is split into ``J`` intervals and the probability for
+    each interval is computed as the difference between the value of the triangular
+    distribution function for the interval boundaries. The probability for the first
+    interval is computed as the value of the triangular distribution function for the
+    first interval boundary.
+
+    The triangular distribution function is denoted as :math:`\\text{f}(x, a, b, c)`.
+    The parameters :math:`a`, :math:`b` and :math:`c` for class :math:`j` are computed
+    as described in :footcite:t:`vargas2023gentri`.
 
     Parameters
     ----------
-    n : int
+    J : int
         Number of classes.
     alphas : np.ndarray
         Array of alphas.
     verbose : int, optional
         Verbosity level, by default 0.
+
+    Raises
+    ------
+    ValueError
+        If ``J`` is not a positive integer greater than 1.
+        If ``alphas`` is not a numpy array of shape :math:`(2J,)`.
+
+    Returns
+    -------
+    probs : 2d array-like of shape (J, J)
+        Matrix of probabilities where each row represents the true class
+        and each column the probability for class j.
+
+    Example
+    -------
+    >>> from dlordinal.distributions import get_general_triangular_probabilities
+    >>> get_general_triangular_probabilities(
+    ...     5,
+    ...     [0, 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.1, 0]
+    ... )
+    array([[0.95, 0.05, 0.  , 0.  , 0.  ],
+           [0.05, 0.9 , 0.05, 0.  , 0.  ],
+           [0.  , 0.05, 0.9 , 0.05, 0.  ],
+           [0.  , 0.  , 0.05, 0.85, 0.1 ],
+           [0.  , 0.  , 0.  , 0.  , 1.  ]], dtype=float32)
     """
 
-    intervals = get_intervals(n)
+    if J < 2:
+        raise ValueError(f"J must be greater than 1, but got {J}")
+
+    if isinstance(alphas, list):
+        alphas = np.array(alphas)
+
+    if not isinstance(alphas, np.ndarray) or alphas.shape != (2 * J,):
+        raise ValueError(
+            "alphas must be a numpy array or list of shape (2 * n,),"
+            " but got {alphas.shape}"
+        )
+
+    intervals = get_intervals(J)
     probs = []
 
     # Compute probability for each interval (class) using the distribution function.
-    params = get_general_triangular_params(n, alphas, verbose)
+    params = get_general_triangular_params(J, alphas, verbose)
 
     for param in params:
         a = param["a"]
diff --git a/dlordinal/distributions/poisson_distribution.py b/dlordinal/distributions/poisson_distribution.py
index 1efb1dc..bc4f54e 100644
--- a/dlordinal/distributions/poisson_distribution.py
+++ b/dlordinal/distributions/poisson_distribution.py
@@ -3,24 +3,55 @@
 from scipy.stats import poisson
 
 
-def get_poisson_probabilities(n):
-    """Get probabilities from poisson distribution for n classes.
+def get_poisson_probabilities(J):
+    """Get probabilities from poisson distribution for ``J`` classes or splits using the
+    methodology described in :footcite:t:`liu2020unimodal`.
+    The :math:`[0,1]` interval is split into ``J`` intervals and the probability for
+    each interval is computed as the difference between the value of the poisson
+    probability function for the interval boundaries. The probability for the first
+    interval is computed as the value of the poisson probability function for the first
+    interval boundary. Then, a softmax function is applied to each row of the resulting
+    matrix to obtain valid probabilities.
+
+    The poisson probability function is denoted as :math:`\\text{p}(k, \\lambda)`
+    where :math:`k` is given by the order of the class for which the probability is
+    computed, and :math:`\\lambda` is given by :math:`k` where :math:`k` is the order
+    of the target class.
 
     Parameters
     ----------
-    n : int
-        Number of classes.
+    J : int
+        Number of classes or splits.
+
+    Raises
+    ------
+    ValueError
+        If ``J`` is not a positive integer greater than 1.
 
     Returns
     -------
-    probs : 2d array-like
+    probs : 2d array-like of shape (J, J)
         Matrix of probabilities where each row represents the true class
-        and each column the probability for class n.
+        and each column the probability for class j.
+
+    Example
+    -------
+    >>> from dlordinal.distributions import get_poisson_probabilities
+    >>> get_poisson_probabilities(5)
+    array([[0.23414552, 0.23414552, 0.19480578, 0.17232403, 0.16457916],
+        [0.18896888, 0.21635436, 0.21635436, 0.19768881, 0.18063359],
+        [0.17822335, 0.19688341, 0.21214973, 0.21214973, 0.20059378],
+        [0.17919028, 0.18931175, 0.20370191, 0.21389803, 0.21389803],
+        [0.18400408, 0.18903075, 0.19882883, 0.21031236, 0.21782399]])
+
     """
 
+    if J < 2 or not isinstance(J, int):
+        raise ValueError(f"{J=} must be a positive integer greater than 1")
+
     probs = []
 
-    for true_class in range(1, n + 1):
-        probs.append(poisson.pmf(np.arange(0, n), true_class))
+    for true_class in range(1, J + 1):
+        probs.append(poisson.pmf(np.arange(0, J), true_class))
 
     return softmax(np.array(probs), axis=1)
diff --git a/dlordinal/distributions/tests/test_beta_distribution.py b/dlordinal/distributions/tests/test_beta_distribution.py
index 61093aa..0224164 100644
--- a/dlordinal/distributions/tests/test_beta_distribution.py
+++ b/dlordinal/distributions/tests/test_beta_distribution.py
@@ -1,7 +1,8 @@
 import numpy as np
 import pytest
 
-from ..beta_distribution import beta_dist, beta_func, get_beta_probabilities
+from dlordinal.distributions import get_beta_softlabels
+from dlordinal.distributions.beta_distribution import beta_dist, beta_func
 
 
 def test_beta_inc():
@@ -57,7 +58,7 @@ def test_beta_probabilities():
     p = 2.0
     q = 3.0
     a = 1.0
-    result = get_beta_probabilities(n, p, q, a)
+    result = get_beta_softlabels(n, p, q, a)
     expected_result = [
         0.1808000009216,
         0.34399999942400017,
@@ -74,7 +75,7 @@ def test_beta_probabilities():
     p = 1.5
     q = 2.5
     a = 2.0
-    result = get_beta_probabilities(n, p, q, a)
+    result = get_beta_softlabels(n, p, q, a)
     expected_result = [
         0.05010107325697135,
         0.283232260076362,
diff --git a/dlordinal/distributions/tests/test_binomial_distribution.py b/dlordinal/distributions/tests/test_binomial_distribution.py
index e2610b0..7802221 100644
--- a/dlordinal/distributions/tests/test_binomial_distribution.py
+++ b/dlordinal/distributions/tests/test_binomial_distribution.py
@@ -1,13 +1,13 @@
 import numpy as np
 import pytest
 
-from ..binomial_distribution import get_binomial_probabilities
+from dlordinal.distributions import get_binomial_softlabels
 
 
 def test_get_binomial_probabilities():
     # Case 1: n = 5
     n = 5
-    result = get_binomial_probabilities(n)
+    result = get_binomial_softlabels(n)
     expected_result = np.array(
         [
             [6.561e-01, 2.916e-01, 4.860e-02, 3.600e-03, 1.000e-04],
diff --git a/dlordinal/distributions/tests/test_exponential_distribution.py b/dlordinal/distributions/tests/test_exponential_distribution.py
index 4a9953f..76f2c8d 100644
--- a/dlordinal/distributions/tests/test_exponential_distribution.py
+++ b/dlordinal/distributions/tests/test_exponential_distribution.py
@@ -1,14 +1,14 @@
 import numpy as np
 import pytest
 
-from ..exponential_distribution import get_exponential_probabilities
+from dlordinal.distributions import get_exponential_softlabels
 
 
 def test_get_exponential_probabilities():
     n = 5
     p = 1.0
     tau = 1.0
-    result = get_exponential_probabilities(n, p, tau)
+    result = get_exponential_softlabels(n, p, tau)
     expected_result = np.array(
         [
             [0.63640865, 0.23412166, 0.08612854, 0.03168492, 0.01165623],
@@ -35,7 +35,7 @@ def test_exponential_probabilities():
     n = 4
     p = 2.0
     tau = 1.0
-    result = get_exponential_probabilities(n, p, tau)
+    result = get_exponential_softlabels(n, p, tau)
     expected_result = np.array(
         [
             [7.21334965e-01, 2.65364304e-01, 1.32117107e-02, 8.90198068e-05],
diff --git a/dlordinal/distributions/tests/test_general_triangular_distribution.py b/dlordinal/distributions/tests/test_general_triangular_distribution.py
index 2b3b8ac..390e1ff 100644
--- a/dlordinal/distributions/tests/test_general_triangular_distribution.py
+++ b/dlordinal/distributions/tests/test_general_triangular_distribution.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from ..general_triangular_distribution import get_general_triangular_params
+from dlordinal.distributions import get_general_triangular_params
 
 
 def test_get_general_triangular_params():
diff --git a/dlordinal/distributions/tests/test_poisson_distribution.py b/dlordinal/distributions/tests/test_poisson_distribution.py
index 6644113..2d4df26 100644
--- a/dlordinal/distributions/tests/test_poisson_distribution.py
+++ b/dlordinal/distributions/tests/test_poisson_distribution.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from ..poisson_distribution import get_poisson_probabilities
+from dlordinal.distributions import get_poisson_probabilities
 
 
 def test_get_poisson_probabilities():
diff --git a/dlordinal/distributions/tests/test_triangular_distribution.py b/dlordinal/distributions/tests/test_triangular_distribution.py
index 1391b7d..511a4f9 100644
--- a/dlordinal/distributions/tests/test_triangular_distribution.py
+++ b/dlordinal/distributions/tests/test_triangular_distribution.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from ..triangular_distribution import get_triangular_probabilities
+from dlordinal.distributions import get_triangular_softlabels
 
 
 def test_get_triangular_probabilities():
@@ -10,7 +10,7 @@ def test_get_triangular_probabilities():
     alpha2 = 0.01
     verbose = 0
 
-    result = get_triangular_probabilities(n, alpha2, verbose)
+    result = get_triangular_softlabels(n, alpha2, verbose)
 
     expected_result = [
         [0.98845494, 0.01154505, 0.0, 0.0, 0.0],
@@ -39,7 +39,7 @@ def test_get_triangular_probabilities():
     alpha2 = 0.01
     verbose = 0
 
-    result = get_triangular_probabilities(n, alpha2, verbose)
+    result = get_triangular_softlabels(n, alpha2, verbose)
 
     expected_result = [
         [0.98845494, 0.01154505, 0.0, 0.0, 0.0, 0.0, 0.0],
@@ -72,7 +72,7 @@ def test_get_triangular_probabilities_verbose():
     alpha2 = 0.01
     verbose = 4
 
-    result = get_triangular_probabilities(n, alpha2, verbose)
+    result = get_triangular_softlabels(n, alpha2, verbose)
 
     expected_result = [
         [
diff --git a/dlordinal/distributions/tests/test_utils.py b/dlordinal/distributions/tests/test_utils.py
index d246d99..2bd1c8c 100644
--- a/dlordinal/distributions/tests/test_utils.py
+++ b/dlordinal/distributions/tests/test_utils.py
@@ -1,4 +1,4 @@
-from ..utils import get_intervals, triangular_cdf
+from dlordinal.distributions.utils import get_intervals, triangular_cdf
 
 
 def test_get_intervals():
diff --git a/dlordinal/distributions/triangular_distribution.py b/dlordinal/distributions/triangular_distribution.py
index 4988123..4deb6ef 100644
--- a/dlordinal/distributions/triangular_distribution.py
+++ b/dlordinal/distributions/triangular_distribution.py
@@ -5,22 +5,105 @@
 from .utils import get_intervals, triangular_cdf
 
 
-def get_triangular_probabilities(n: int, alpha2: float = 0.01, verbose: int = 0):
+def get_triangular_softlabels(J: int, alpha2: float = 0.01, verbose: int = 0):
     """
-    Get the probabilities for the triangular distribution.
+    Get probabilities from triangular distributions for ``J`` classes or splits using
+    the approach described in :footcite:t:`vargas2023softlabelling`.
+    The :math:`[0,1]` interval is split into ``J`` intervals and the probability for
+    each interval is computed as the difference between the value of the triangular
+    distribution function for the interval boundaries. The probability for the first
+    interval is computed as the value of the triangular distribution function for the
+    first interval boundary.
+
+    The triangular distribution function is denoted as :math:`\\text{p}(x, a, b, c)`
+    where :math:`a`, :math:`b` and :math:`c` are the parameters of the distribution,
+    and are determined by the number of classes :math:`J` and the value of the
+    :math:`\\alpha_2` parameter. The value of :math:`\\alpha_2` represents the
+    probability that is assigned to the adjacent classes of the target class. The
+    parameters :math:`a`, :math:`b` and :math:`c` for class :math:`j` are computed
+    as follows:
+
+    .. math::
+        a_j = \\begin{cases}
+            0 & \\text{if } j = 1 \\\\
+            \\frac{2j - 2 - 4j\\alpha_2 + 2\\alpha_2 \\pm
+              \\sqrt{2\\alpha_2}}{2n(1 - 2\\alpha_2)} & \\text{if } 1 < j < J\\\\
+            1 + \\frac{1}{\\pm J (\\sqrt{\\alpha_3} - 1)} & \\text{if } j = J
+        \\end{cases}
+
+    .. math::
+        b_j = \\begin{cases}
+            \\frac{1}{(1 - \\sqrt{\\alpha_1})n} & \\text{if } j = 1 \\\\
+            \\frac{2j - 4j\\alpha_2 + 2\\alpha_2 \\pm
+             \\sqrt{2\\alpha_2}}{2n(1 - 2\\alpha_2)} & \\text{if } 1 < j < J\\\\
+             1 & \\text{if } j = J
+        \\end{cases}
+
+    .. math::
+        c_j = \\begin{cases}
+            0 & \\text{if } j = 1 \\\\
+            \\frac{a + b}{2} & \\text{if } 1 < j  < J\\\\
+            1 & \\text{if } j = J
+        \\end{cases}
+
+    The value of :math:`\\alpha_1`, that represents the error for the first class,
+    is computed as follows:
+
+    .. math::
+        \\alpha_1 = \\left(\\frac{1 - \\sqrt{1 - 4(1 - 2\\alpha_2)(2\\alpha_2 -
+          \\sqrt{2\\alpha_2})}}{2}\\right)^2
+
+    The value of :math:`\\alpha_3`, that represents the error for the last class,
+    is computed as follows:
+
+    .. math::
+        \\alpha_3 = \\left(\\frac{1 -
+        \\sqrt{1 - 4\\left(\\frac{J - 1}{J}\\right)^2(1 - 2\\alpha_2)(\\sqrt{2\\alpha_2}
+          (-1 + \\sqrt{2\\alpha_2}))}}{2}\\right)^2
+
+
+    The value of :math:`\\alpha_2` is given by the user.
 
     Parameters
     ----------
-    n : int
-        Number of classes.
-    alpha2 : float, optional
-        Alpha2 value, by default 0.01.
-    verbose : int, optional
-        Verbosity level, by default 0.
+    J : int
+        Number of classes or splits (:math:`J`).
+    alpha2 : float, optional, default=0.01
+        Value of the :math:`\\alpha_2` parameter.
+    verbose : int, optional, default=0
+        Verbosity level.
+
+    Raises
+    ------
+    ValueError
+        If ``J`` is not a positive integer greater than 1.
+        If ``alpha2`` is not a float between 0 and 1.
+
+    Returns
+    -------
+    probs : 2d array-like of shape (J, J)
+        Matrix of probabilities where each row represents the true class
+        and each column the probability for class j.
+
+    Example
+    -------
+    >>> from dlordinal.distributions import get_triangular_probabilities
+    >>> get_triangular_probabilities(5)
+    array([[0.98845494, 0.01154505, 0.        , 0.        , 0.        ],
+           [0.01      , 0.98      , 0.01      , 0.        , 0.        ],
+           [0.        , 0.01      , 0.98      , 0.01      , 0.        ],
+           [0.        , 0.        , 0.01      , 0.98      , 0.01      ],
+           [0.        , 0.        , 0.        , 0.00505524, 0.99494475]])
     """
 
+    if J < 2 or not isinstance(J, int):
+        raise ValueError(f"{J=} must be a positive integer greater than 1")
+
+    if alpha2 < 0 or alpha2 > 1:
+        raise ValueError(f"{alpha2=} must be a float between 0 and 1")
+
     if verbose >= 1:
-        print(f"Computing triangular probabilities for {n=} and {alpha2=}...")
+        print(f"Computing triangular probabilities for {J=} and {alpha2=}...")
 
     def compute_alpha1(alpha2):
         c_minus = (1 - 2 * alpha2) * (2 * alpha2 - math.sqrt(2 * alpha2))
@@ -29,7 +112,7 @@ def compute_alpha1(alpha2):
 
     def compute_alpha3(alpha2):
         c1 = (
-            pow((n - 1) / n, 2)
+            pow((J - 1) / J, 2)
             * (1 - 2 * alpha2)
             * (math.sqrt(2 * alpha2) * (-1 + math.sqrt(2 * alpha2)))
         )
@@ -42,18 +125,15 @@ def compute_alpha3(alpha2):
     if verbose >= 1:
         print(f"{alpha1=}, {alpha2=}, {alpha3=}")
 
-    def b1(n):
-        return 1.0 / ((1.0 - math.sqrt(alpha1)) * n)
+    def b1(J):
+        return 1.0 / ((1.0 - math.sqrt(alpha1)) * J)
 
-    def m1(n):
-        return math.sqrt(alpha1) / ((1 - math.sqrt(alpha1)) * n)
-
-    def aj(n, j):
+    def aj(J, j):
         num1 = 2.0 * j - 2 - 4 * j * alpha2 + 2 * alpha2
         num2 = math.sqrt(2 * alpha2)
-        den = 2.0 * n * (1 - 2 * alpha2)
+        den = 2.0 * J * (1 - 2 * alpha2)
 
-        max_value = (j - 1.0) / n
+        max_value = (j - 1.0) / J
 
         # +-
         return (
@@ -62,12 +142,12 @@ def aj(n, j):
             else (num1 - num2) / den
         )
 
-    def bj(n, j):
+    def bj(J, j):
         num1 = 2.0 * j - 4 * j * alpha2 + 2 * alpha2
         num2 = math.sqrt(2 * alpha2)
-        den = 2.0 * n * (1 - 2 * alpha2)
+        den = 2.0 * J * (1 - 2 * alpha2)
 
-        min_value = j / n
+        min_value = j / J
 
         # +-
         return (
@@ -76,42 +156,37 @@ def bj(n, j):
             else (num1 - num2) / den
         )
 
-    def aJ(n):
-        aJ_plus = 1.0 + 1.0 / (n * (math.sqrt(alpha3) - 1.0))
-        aJ_minus = 1.0 + 1.0 / (-n * (math.sqrt(alpha3) - 1.0))
+    def aJ(J):
+        aJ_plus = 1.0 + 1.0 / (J * (math.sqrt(alpha3) - 1.0))
+        aJ_minus = 1.0 + 1.0 / (-J * (math.sqrt(alpha3) - 1.0))
         return aJ_plus if aJ_plus > 0.0 else aJ_minus
 
-    def nJ(n):
-        num = math.sqrt(alpha3)
-        den = n * (1 - alpha3)
-        return num / den
-
     if verbose >= 3:
-        print(f"{b1(n)=}, {m1(n)=}, {aJ(n)=}, {nJ(n)=}, {aj(n, 1)=}, {bj(n,1)=}")
-        for i in range(1, n + 1):
-            print(f"{i=}  {aj(n, i)=}, {bj(n,i)=}")
+        print(f"{b1(J)=}, {aJ(J)=}, {aj(J, 1)=}, {bj(J,1)=}")
+        for i in range(1, J + 1):
+            print(f"{i=}  {aj(J, i)=}, {bj(J,i)=}")
 
-    intervals = get_intervals(n)
+    intervals = get_intervals(J)
     probs = []
 
     # Compute probability for each interval (class) using the distribution function.
-    for j in range(1, n + 1):
+    for j in range(1, J + 1):
         j_probs = []
         if j == 1:
             a = 0.0
-            b = b1(n)
+            b = b1(J)
             c = 0.0
-        elif j == n:
-            a = aJ(n)
+        elif j == J:
+            a = aJ(J)
             b = 1.0
             c = 1.0
         else:
-            a = aj(n, j)
-            b = bj(n, j)
+            a = aj(J, j)
+            b = bj(J, j)
             c = (a + b) / 2.0
 
         if verbose >= 1:
-            print(f"Class: {j}, {a=}, {b=}, {c=}, (j-1)/J={(j-1)/n}, (j/J)={j/n}")
+            print(f"Class: {j}, {a=}, {b=}, {c=}, (j-1)/J={(j-1)/J}, (j/J)={j/J}")
 
         for interval in intervals:
             j_probs.append(
diff --git a/dlordinal/estimator/pytorch_estimator.py b/dlordinal/estimator/pytorch_estimator.py
index 17e6b34..311d40e 100644
--- a/dlordinal/estimator/pytorch_estimator.py
+++ b/dlordinal/estimator/pytorch_estimator.py
@@ -1,5 +1,6 @@
 from typing import Optional, Union
 
+import numpy as np
 import torch
 import torch.nn.functional as F
 from sklearn.base import BaseEstimator
@@ -166,7 +167,7 @@ def predict_proba(self, X: Union[DataLoader, torch.Tensor]):
                 predictions.append(predictions_batch)
 
             # Concatenate predictions
-            predictions = torch.cat(predictions)
+            predictions = np.concatenate(predictions)
             return predictions
 
         # check if X is a torch Tensor
@@ -191,7 +192,7 @@ def _predict_proba(self, X):
             X = X.to(self.device)
             pred = self.model(X)
             probabilities = F.softmax(pred, dim=1)
-            return probabilities
+            return probabilities.cpu().numpy()
 
     def predict(self, X: Union[DataLoader, torch.Tensor]):
         """
@@ -203,4 +204,4 @@ def predict(self, X: Union[DataLoader, torch.Tensor]):
             The data to predict.
         """
         pred = self.predict_proba(X)
-        return torch.argmax(pred, dim=1)
+        return np.argmax(pred, axis=1)
diff --git a/dlordinal/estimator/tests/test_estimator.py b/dlordinal/estimator/tests/test_estimator.py
index 91b4bf2..f73d7ad 100644
--- a/dlordinal/estimator/tests/test_estimator.py
+++ b/dlordinal/estimator/tests/test_estimator.py
@@ -4,7 +4,7 @@
 from torch.utils.data import DataLoader, TensorDataset
 from torchvision import models
 
-from ..pytorch_estimator import PytorchEstimator
+from dlordinal.estimator import PytorchEstimator
 
 
 def test_pytorch_estimator_creation():
@@ -139,7 +139,7 @@ def test_pytorch_estimator_predict():
     assert len(predictions) == 50
 
     # Check that the predictions are values in the range [0, num_classes)
-    assert torch.all(predictions >= 0) and torch.all(predictions < num_classes)
+    assert np.all(predictions >= 0) and np.all(predictions < num_classes)
 
 
 def test_pytorch_estimator_predict_proba_dataloader():
@@ -175,10 +175,10 @@ def test_pytorch_estimator_predict_proba_dataloader():
     assert probabilities.shape == (50, 5)
 
     # Verify that the sum of the probabilities for each example is close to 1.
-    assert torch.allclose(probabilities.sum(dim=1), torch.ones(50), atol=1e-5)
+    assert np.allclose(np.sum(probabilities, axis=1), np.ones(50), atol=1e-5)
 
     # Verify that the probabilities are in the range [0, 1]
-    assert torch.all(probabilities >= 0) and torch.all(probabilities <= 1)
+    assert np.all(probabilities >= 0) and np.all(probabilities <= 1)
 
 
 def test_pytorch_estimator_predict_proba_tensor():
@@ -214,7 +214,7 @@ def test_pytorch_estimator_predict_proba_tensor():
     assert probabilities.shape == (50, 3)
 
     # Check that the sum of the probabilities for each example is close to 1.
-    assert torch.allclose(probabilities.sum(dim=1), torch.ones(50), atol=1e-5)
+    assert np.allclose(np.sum(probabilities, axis=1), np.ones(50), atol=1e-5)
 
     # Check that the probabilities are in the range [0, 1]
-    assert torch.all(probabilities >= 0) and torch.all(probabilities <= 1)
+    assert np.all(probabilities >= 0) and np.all(probabilities <= 1)
diff --git a/dlordinal/estimator/tests/test_types.py b/dlordinal/estimator/tests/test_types.py
index 8b18f7e..4c691f8 100644
--- a/dlordinal/estimator/tests/test_types.py
+++ b/dlordinal/estimator/tests/test_types.py
@@ -2,7 +2,7 @@
 import torch
 from torch.utils.data import DataLoader, TensorDataset
 
-from ..pytorch_estimator import PytorchEstimator
+from dlordinal.estimator import PytorchEstimator
 
 
 @pytest.fixture
diff --git a/dlordinal/layers/tests/test_clm.py b/dlordinal/layers/tests/test_clm.py
index cc07e6e..0266e0d 100644
--- a/dlordinal/layers/tests/test_clm.py
+++ b/dlordinal/layers/tests/test_clm.py
@@ -1,6 +1,6 @@
 import torch
 
-from ..clm import CLM
+from dlordinal.layers import CLM
 
 
 def test_clm_creation():
diff --git a/dlordinal/layers/tests/test_ordinal_fully_connected.py b/dlordinal/layers/tests/test_ordinal_fully_connected.py
index 7144cbc..4036b86 100644
--- a/dlordinal/layers/tests/test_ordinal_fully_connected.py
+++ b/dlordinal/layers/tests/test_ordinal_fully_connected.py
@@ -1,10 +1,10 @@
 import pytest
 import torch
 
-from ..activation_function import activation_function_by_name
-from ..ordinal_fully_connected import (
+from dlordinal.layers import (
     ResNetOrdinalFullyConnected,
     VGGOrdinalFullyConnected,
+    activation_function_by_name,
 )
 
 
diff --git a/dlordinal/layers/tests/test_stick_breaking_layer.py b/dlordinal/layers/tests/test_stick_breaking_layer.py
index 657c3fa..59b4471 100644
--- a/dlordinal/layers/tests/test_stick_breaking_layer.py
+++ b/dlordinal/layers/tests/test_stick_breaking_layer.py
@@ -1,6 +1,6 @@
 import torch
 
-from ..stick_breaking_layer import StickBreakingLayer
+from dlordinal.layers import StickBreakingLayer
 
 
 def test_stick_breaking_layer_creation():
diff --git a/dlordinal/losses/beta_loss.py b/dlordinal/losses/beta_loss.py
index ed78e26..30dfa09 100644
--- a/dlordinal/losses/beta_loss.py
+++ b/dlordinal/losses/beta_loss.py
@@ -3,7 +3,7 @@
 import torch
 from torch import Tensor
 
-from ..distributions import get_beta_probabilities
+from ..distributions import get_beta_softlabels
 from .custom_targets_loss import CustomTargetsCrossEntropyLoss
 
 # Params [a,b] for beta distribution
@@ -177,7 +177,7 @@ def __init__(
         # Precompute class probabilities for each label
         cls_probs = torch.tensor(
             [
-                get_beta_probabilities(
+                get_beta_softlabels(
                     num_classes,
                     self.params[num_classes][i][0],
                     self.params[num_classes][i][1],
diff --git a/dlordinal/losses/binomial_loss.py b/dlordinal/losses/binomial_loss.py
index 0fd88c4..1c1f51a 100644
--- a/dlordinal/losses/binomial_loss.py
+++ b/dlordinal/losses/binomial_loss.py
@@ -3,7 +3,7 @@
 import torch
 from torch import Tensor
 
-from ..distributions import get_binomial_probabilities
+from ..distributions import get_binomial_softlabels
 from .custom_targets_loss import CustomTargetsCrossEntropyLoss
 
 
@@ -59,7 +59,7 @@ def __init__(
         label_smoothing: float = 0.0,
     ):
         # Precompute class probabilities for each label
-        cls_probs = torch.tensor(get_binomial_probabilities(num_classes)).float()
+        cls_probs = torch.tensor(get_binomial_softlabels(num_classes)).float()
 
         super().__init__(
             cls_probs=cls_probs,
diff --git a/dlordinal/losses/exponential_loss.py b/dlordinal/losses/exponential_loss.py
index e2f6b71..e0fecc5 100644
--- a/dlordinal/losses/exponential_loss.py
+++ b/dlordinal/losses/exponential_loss.py
@@ -3,7 +3,7 @@
 import torch
 from torch import Tensor
 
-from ..distributions import get_exponential_probabilities
+from ..distributions import get_exponential_softlabels
 from .custom_targets_loss import CustomTargetsCrossEntropyLoss
 
 
@@ -63,7 +63,7 @@ def __init__(
         label_smoothing: float = 0.0,
     ):
         # Precompute class probabilities for each label
-        cls_probs = torch.tensor(get_exponential_probabilities(num_classes, p)).float()
+        cls_probs = torch.tensor(get_exponential_softlabels(num_classes, p)).float()
 
         super().__init__(
             cls_probs=cls_probs,
diff --git a/dlordinal/losses/general_triangular_loss.py b/dlordinal/losses/general_triangular_loss.py
index 3de6ed6..0ad1c06 100644
--- a/dlordinal/losses/general_triangular_loss.py
+++ b/dlordinal/losses/general_triangular_loss.py
@@ -4,7 +4,7 @@
 import torch
 from torch import Tensor
 
-from ..distributions import get_general_triangular_probabilities
+from ..distributions import get_general_triangular_softlabels
 from .custom_targets_loss import CustomTargetsCrossEntropyLoss
 
 
@@ -63,7 +63,7 @@ def __init__(
         label_smoothing: float = 0.0,
     ):
         # Precompute class probabilities for each label
-        r = get_general_triangular_probabilities(num_classes, alphas, verbose=0)
+        r = get_general_triangular_softlabels(num_classes, alphas, verbose=0)
         cls_probs = torch.tensor(r)
 
         super().__init__(
diff --git a/dlordinal/losses/tests/test_beta_loss.py b/dlordinal/losses/tests/test_beta_loss.py
index 4fda664..b48b2a5 100644
--- a/dlordinal/losses/tests/test_beta_loss.py
+++ b/dlordinal/losses/tests/test_beta_loss.py
@@ -1,7 +1,7 @@
 import pytest
 import torch
 
-from ..beta_loss import BetaCrossEntropyLoss
+from dlordinal.losses import BetaCrossEntropyLoss
 
 
 def test_beta_loss_creation():
diff --git a/dlordinal/losses/tests/test_binomial_loss.py b/dlordinal/losses/tests/test_binomial_loss.py
index 1605016..61dec39 100644
--- a/dlordinal/losses/tests/test_binomial_loss.py
+++ b/dlordinal/losses/tests/test_binomial_loss.py
@@ -1,7 +1,7 @@
 import pytest
 import torch
 
-from ..binomial_loss import BinomialCrossEntropyLoss
+from dlordinal.losses import BinomialCrossEntropyLoss
 
 
 def test_binomial_loss_creation():
diff --git a/dlordinal/losses/tests/test_custom_loss.py b/dlordinal/losses/tests/test_custom_loss.py
index 730863f..b4802e6 100644
--- a/dlordinal/losses/tests/test_custom_loss.py
+++ b/dlordinal/losses/tests/test_custom_loss.py
@@ -1,6 +1,6 @@
 import torch
 
-from ..custom_targets_loss import CustomTargetsCrossEntropyLoss
+from dlordinal.losses import CustomTargetsCrossEntropyLoss
 
 
 # Auxiliar function to get a test class
diff --git a/dlordinal/losses/tests/test_exponential_loss.py b/dlordinal/losses/tests/test_exponential_loss.py
index b41972d..b5536bf 100644
--- a/dlordinal/losses/tests/test_exponential_loss.py
+++ b/dlordinal/losses/tests/test_exponential_loss.py
@@ -1,7 +1,7 @@
 import pytest
 import torch
 
-from ..exponential_loss import ExponentialRegularisedCrossEntropyLoss
+from dlordinal.losses import ExponentialRegularisedCrossEntropyLoss
 
 
 def test_exponential_loss_creation():
diff --git a/dlordinal/losses/tests/test_general_triangular_loss.py b/dlordinal/losses/tests/test_general_triangular_loss.py
index fe5d50e..0227bd9 100644
--- a/dlordinal/losses/tests/test_general_triangular_loss.py
+++ b/dlordinal/losses/tests/test_general_triangular_loss.py
@@ -2,7 +2,7 @@
 import pytest
 import torch
 
-from ..general_triangular_loss import GeneralTriangularCrossEntropyLoss
+from dlordinal.losses import GeneralTriangularCrossEntropyLoss
 
 
 def test_general_triangular_loss_creation():
diff --git a/dlordinal/losses/tests/test_mceloss.py b/dlordinal/losses/tests/test_mceloss.py
index f05fe78..c9283d2 100644
--- a/dlordinal/losses/tests/test_mceloss.py
+++ b/dlordinal/losses/tests/test_mceloss.py
@@ -1,7 +1,7 @@
 import pytest
 import torch
 
-from ..mceloss import MCELoss
+from dlordinal.losses import MCELoss
 
 
 def test_mceloss_creation():
diff --git a/dlordinal/losses/tests/test_mcewkloss.py b/dlordinal/losses/tests/test_mcewkloss.py
index 5fff804..16daef5 100644
--- a/dlordinal/losses/tests/test_mcewkloss.py
+++ b/dlordinal/losses/tests/test_mcewkloss.py
@@ -1,9 +1,7 @@
 import pytest
 import torch
 
-from ..mceloss import MCELoss
-from ..mcewkloss import MCEAndWKLoss
-from ..wkloss import WKLoss
+from dlordinal.losses import MCEAndWKLoss, MCELoss, WKLoss
 
 
 def test_mcewkloss_creation():
diff --git a/dlordinal/losses/tests/test_ordinal_ecoc_loss.py b/dlordinal/losses/tests/test_ordinal_ecoc_loss.py
index d8c85dd..b41f632 100644
--- a/dlordinal/losses/tests/test_ordinal_ecoc_loss.py
+++ b/dlordinal/losses/tests/test_ordinal_ecoc_loss.py
@@ -1,7 +1,7 @@
 import torch
 from torch import cuda
 
-from ..ordinal_ecoc_distance_loss import OrdinalECOCDistanceLoss
+from dlordinal.losses import OrdinalECOCDistanceLoss
 
 
 def test_ordinal_ecoc_distance_loss_creation():
diff --git a/dlordinal/losses/tests/test_poisson_loss.py b/dlordinal/losses/tests/test_poisson_loss.py
index f03c83c..a7d87de 100644
--- a/dlordinal/losses/tests/test_poisson_loss.py
+++ b/dlordinal/losses/tests/test_poisson_loss.py
@@ -1,7 +1,7 @@
 import pytest
 import torch
 
-from ..poisson_loss import PoissonCrossEntropyLoss
+from dlordinal.losses import PoissonCrossEntropyLoss
 
 
 def test_poisson_loss_creation():
diff --git a/dlordinal/losses/tests/test_triangular_loss.py b/dlordinal/losses/tests/test_triangular_loss.py
index 6ebc4a0..0f0af72 100644
--- a/dlordinal/losses/tests/test_triangular_loss.py
+++ b/dlordinal/losses/tests/test_triangular_loss.py
@@ -1,7 +1,7 @@
 import pytest
 import torch
 
-from ..triangular_loss import TriangularCrossEntropyLoss
+from dlordinal.losses import TriangularCrossEntropyLoss
 
 
 def test_triangular_loss_creation():
diff --git a/dlordinal/losses/tests/test_wkloss.py b/dlordinal/losses/tests/test_wkloss.py
index 9b1055f..89411a8 100644
--- a/dlordinal/losses/tests/test_wkloss.py
+++ b/dlordinal/losses/tests/test_wkloss.py
@@ -1,7 +1,7 @@
 import pytest
 import torch
 
-from ..wkloss import WKLoss
+from dlordinal.losses import WKLoss
 
 
 def test_wkloss_creation():
diff --git a/dlordinal/losses/triangular_loss.py b/dlordinal/losses/triangular_loss.py
index b25e873..2e5c2b7 100644
--- a/dlordinal/losses/triangular_loss.py
+++ b/dlordinal/losses/triangular_loss.py
@@ -3,7 +3,7 @@
 import torch
 from torch import Tensor
 
-from ..distributions import get_triangular_probabilities
+from ..distributions import get_triangular_softlabels
 from .custom_targets_loss import CustomTargetsCrossEntropyLoss
 
 
@@ -62,7 +62,7 @@ def __init__(
         label_smoothing: float = 0.0,
     ):
         # Precompute class probabilities for each label
-        cls_probs = torch.tensor(get_triangular_probabilities(num_classes, alpha2))
+        cls_probs = torch.tensor(get_triangular_softlabels(num_classes, alpha2))
         super().__init__(
             cls_probs=cls_probs,
             eta=eta,
diff --git a/dlordinal/metrics/tests/test_metrics.py b/dlordinal/metrics/tests/test_metrics.py
index 8b17ead..f8914f3 100644
--- a/dlordinal/metrics/tests/test_metrics.py
+++ b/dlordinal/metrics/tests/test_metrics.py
@@ -5,7 +5,7 @@
 import pytest
 from sklearn.metrics import recall_score
 
-from ..metrics import (
+from dlordinal.metrics import (
     accuracy_off1,
     gmsec,
     minimum_sensitivity,
diff --git a/docs/distributions.rst b/docs/distributions.rst
index d49ffb6..fd0a674 100644
--- a/docs/distributions.rst
+++ b/docs/distributions.rst
@@ -5,3 +5,5 @@ Probability distributions
 
 .. automodule:: dlordinal.distributions
     :members:
+
+.. footbibliography::
diff --git a/tutorials/ecoc_tutorial.ipynb b/tutorials/ecoc_tutorial.ipynb
index 9555460..fcf1ad1 100644
--- a/tutorials/ecoc_tutorial.ipynb
+++ b/tutorials/ecoc_tutorial.ipynb
@@ -73,31 +73,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Downloading http://yanweifu.github.io/FG_NET_data/FGNET.zip to datasets/fgnet/fgnet.zip\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "100%|██████████| 46221653/46221653 [00:01<00:00, 40386000.42it/s]\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Extracting datasets/fgnet/fgnet.zip to datasets/fgnet\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "fgnet = FGNet(root=\"./datasets/fgnet\", download=True, process_data=True)\n",
     "\n",
diff --git a/tutorials/estimator_tutorial.ipynb b/tutorials/estimator_tutorial.ipynb
index b3cfa73..e3f4de3 100644
--- a/tutorials/estimator_tutorial.ipynb
+++ b/tutorials/estimator_tutorial.ipynb
@@ -142,7 +142,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = models.resnet18(pretrained=True)\n",
+    "model = models.resnet18(weights='IMAGENET1K_V1')\n",
     "model.fc = nn.Linear(model.fc.in_features, num_classes)\n",
     "model = model.to(device)\n",
     "\n",
@@ -282,7 +282,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.7"
+   "version": "3.8.18"
   }
  },
  "nbformat": 4,
diff --git a/tutorials/losses_tutorial.ipynb b/tutorials/losses_tutorial.ipynb
index c9358d8..2130239 100644
--- a/tutorials/losses_tutorial.ipynb
+++ b/tutorials/losses_tutorial.ipynb
@@ -178,7 +178,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = models.resnet18(pretrained=True)\n",
+    "model = models.resnet18(weights='IMAGENET1K_V1')\n",
     "model.fc = nn.Linear(model.fc.in_features, num_classes)\n",
     "model = model.to(device)\n",
     "\n",
diff --git a/tutorials/stick_breaking_tutorial.ipynb b/tutorials/stick_breaking_tutorial.ipynb
index b1a05b9..ce05522 100644
--- a/tutorials/stick_breaking_tutorial.ipynb
+++ b/tutorials/stick_breaking_tutorial.ipynb
@@ -171,7 +171,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = models.resnet18(pretrained=True)\n",
+    "model = models.resnet18(weights='IMAGENET1K_V1')\n",
     "model.fc = StickBreakingLayer(model.fc.in_features, num_classes)\n",
     "model = model.to(device)\n",
     "\n",