diff --git a/.gitignore b/.gitignore
index ef3adce17..cabbb8925 100644
--- a/.gitignore
+++ b/.gitignore
@@ -248,3 +248,4 @@ logo/
 /python_model_third_party.py
 /third_party_script.py
 /docs/source/auto_examples/
+scratch/
\ No newline at end of file
diff --git a/src/UQpy/distributions/collection/Beta.py b/src/UQpy/distributions/collection/Beta.py
index f2b5804c0..92f4d0a70 100644
--- a/src/UQpy/distributions/collection/Beta.py
+++ b/src/UQpy/distributions/collection/Beta.py
@@ -1,5 +1,4 @@
 from typing import Union
-
 import scipy.stats as stats
 from UQpy.distributions.baseclass import DistributionContinuous1D
 from beartype import beartype
@@ -30,3 +29,11 @@ def __init__(
             ordered_parameters=("a", "b", "loc", "scale"),
         )
         self._construct_from_scipy(scipy_name=stats.beta)
+
+    def __repr__(self):
+        s = "{a}, {b}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "Beta(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Binomial.py b/src/UQpy/distributions/collection/Binomial.py
index af402bcaa..5a4350ecc 100644
--- a/src/UQpy/distributions/collection/Binomial.py
+++ b/src/UQpy/distributions/collection/Binomial.py
@@ -22,3 +22,9 @@ def __init__(
         """
         super().__init__(n=n, p=p, loc=loc, ordered_parameters=("n", "p", "loc"))
         self._construct_from_scipy(scipy_name=stats.binom)
+
+    def __repr__(self):
+        s = "{n}, {p}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        return "Binomial(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Cauchy.py b/src/UQpy/distributions/collection/Cauchy.py
index 51ce9b24d..ccc974fe2 100644
--- a/src/UQpy/distributions/collection/Cauchy.py
+++ b/src/UQpy/distributions/collection/Cauchy.py
@@ -19,3 +19,12 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.cauchy)
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append("loc={loc}")
+        if self.parameters["scale"] != 1.0:
+            s.append("scale={scale}")
+        s = ", ".join(s)
+        return "Cauchy(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/ChiSquare.py b/src/UQpy/distributions/collection/ChiSquare.py
index dea8bea9d..ccb42eb84 100644
--- a/src/UQpy/distributions/collection/ChiSquare.py
+++ b/src/UQpy/distributions/collection/ChiSquare.py
@@ -24,3 +24,11 @@ def __init__(
             df=df, loc=loc, scale=scale, ordered_parameters=("df", "loc", "scale")
         )
         self._construct_from_scipy(scipy_name=stats.chi2)
+
+    def __repr__(self):
+        s = "{df}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "ChiSquare(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Exponential.py b/src/UQpy/distributions/collection/Exponential.py
index 7376bfba9..1083d046f 100644
--- a/src/UQpy/distributions/collection/Exponential.py
+++ b/src/UQpy/distributions/collection/Exponential.py
@@ -19,3 +19,12 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.expon)
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append("loc={loc}")
+        if self.parameters["scale"] != 1.0:
+            s.append("scale={scale}")
+        s = ", ".join(s)
+        return "Exponential(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Gamma.py b/src/UQpy/distributions/collection/Gamma.py
index 0e5db0744..837c9852f 100644
--- a/src/UQpy/distributions/collection/Gamma.py
+++ b/src/UQpy/distributions/collection/Gamma.py
@@ -25,3 +25,11 @@ def __init__(
             a=a, loc=loc, scale=scale, ordered_parameters=("a", "loc", "scale")
         )
         self._construct_from_scipy(scipy_name=stats.gamma)
+
+    def __repr__(self):
+        s = "{a}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "Gamma(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/GeneralizedExtreme.py b/src/UQpy/distributions/collection/GeneralizedExtreme.py
index 0982d8e76..2f73de944 100644
--- a/src/UQpy/distributions/collection/GeneralizedExtreme.py
+++ b/src/UQpy/distributions/collection/GeneralizedExtreme.py
@@ -25,3 +25,11 @@ def __init__(
             c=c, loc=loc, scale=scale, ordered_parameters=("c", "loc", "scale")
         )
         self._construct_from_scipy(scipy_name=stats.genextreme)
+
+    def __repr__(self):
+        s = "{c}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "GeneralizedExtreme(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/InverseGaussian.py b/src/UQpy/distributions/collection/InverseGaussian.py
index d95315c17..4ead87a0d 100644
--- a/src/UQpy/distributions/collection/InverseGaussian.py
+++ b/src/UQpy/distributions/collection/InverseGaussian.py
@@ -25,3 +25,11 @@ def __init__(
             mu=mu, loc=loc, scale=scale, ordered_parameters=("mu", "loc", "scale")
         )
         self._construct_from_scipy(scipy_name=stats.invgauss)
+
+    def __repr__(self):
+        s = "{mu}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "InverseGauss(" + s.format(**self.parameters) + ")"
\ No newline at end of file
diff --git a/src/UQpy/distributions/collection/JointCopula.py b/src/UQpy/distributions/collection/JointCopula.py
index ce4cf04d0..0140bdfd8 100644
--- a/src/UQpy/distributions/collection/JointCopula.py
+++ b/src/UQpy/distributions/collection/JointCopula.py
@@ -165,3 +165,6 @@ def update_parameters(self, **kwargs: dict):
                 self.copula.parameters[key] = value
             else:
                 self.marginals[int(index)].parameters[key] = value
+
+    def __repr__(self):
+        return f"JointCopula({self.marginals}, {self.copula})"
diff --git a/src/UQpy/distributions/collection/JointIndependent.py b/src/UQpy/distributions/collection/JointIndependent.py
index 8614ec737..ab8e9be7c 100644
--- a/src/UQpy/distributions/collection/JointIndependent.py
+++ b/src/UQpy/distributions/collection/JointIndependent.py
@@ -14,8 +14,8 @@
 class JointIndependent(DistributionND):
     @beartype
     def __init__(
-            self,
-            marginals: Union[list[DistributionContinuous1D], list[DistributionDiscrete1D]],
+        self,
+        marginals: Union[list[DistributionContinuous1D], list[DistributionDiscrete1D]],
     ):
         """
         :param marginals: list of distribution objects that define the marginals.
@@ -24,12 +24,20 @@ def __init__(
         self.ordered_parameters = []
         for i, m in enumerate(marginals):
             self.ordered_parameters.extend(
-                [key + "_" + str(i) for key in m.ordered_parameters])
+                [key + "_" + str(i) for key in m.ordered_parameters]
+            )
 
         # Check and save the marginals
-        if not (isinstance(marginals, list)
-                and all(isinstance(d, (DistributionContinuous1D, DistributionDiscrete1D)) for d in marginals)):
-            raise ValueError("Input marginals must be a list of Distribution1d objects.")
+        if not (
+            isinstance(marginals, list)
+            and all(
+                isinstance(d, (DistributionContinuous1D, DistributionDiscrete1D))
+                for d in marginals
+            )
+        ):
+            raise ValueError(
+                "Input marginals must be a list of Distribution1d objects."
+            )
         self.marginals = marginals
 
         # If all marginals have a method, the joint has it to
@@ -70,6 +78,7 @@ def joint_log_pdf(dist, x):
                 self.log_pmf = MethodType(joint_log_pdf, self)
 
         if all(hasattr(m, "cdf") for m in self.marginals):
+
             def joint_cdf(dist, x):
                 x = dist.check_x_dimension(x)
                 # Compute cdf of independent marginals
@@ -107,8 +116,8 @@ def joint_fit(dist, data):
                 mle_all = {}
                 for ind_m, marg in enumerate(dist.marginals):
                     if any(
-                            param_value is None
-                            for param_value in marg.get_parameters().values()
+                        param_value is None
+                        for param_value in marg.get_parameters().values()
                     ):
                         mle_i = marg.fit(data[:, ind_m])
                     else:
@@ -125,8 +134,15 @@ def joint_fit(dist, data):
             def joint_moments(dist, moments2return="mvsk"):
                 # Go through all marginals
                 if len(moments2return) == 1:
-                    return np.array([marg.moments(moments2return=moments2return) for marg in dist.marginals])
-                moments_ = [np.empty((len(dist.marginals),)) for _ in range(len(moments2return))]
+                    return np.array(
+                        [
+                            marg.moments(moments2return=moments2return)
+                            for marg in dist.marginals
+                        ]
+                    )
+                moments_ = [
+                    np.empty((len(dist.marginals),)) for _ in range(len(moments2return))
+                ]
                 for ind_m, marg in enumerate(dist.marginals):
                     moments_i = marg.moments(moments2return=moments2return)
                     for j in range(len(moments2return)):
@@ -174,3 +190,6 @@ def update_parameters(self, **kwargs: dict):
             key_split = key_indexed.split("_")
             key, index = "_".join(key_split[:-1]), int(key_split[-1])
             self.marginals[index].parameters[key] = value
+
+    def __repr__(self):
+        return f"JointIndependent({self.marginals})"
diff --git a/src/UQpy/distributions/collection/Laplace.py b/src/UQpy/distributions/collection/Laplace.py
index 6c1905cea..0de2d5d1f 100644
--- a/src/UQpy/distributions/collection/Laplace.py
+++ b/src/UQpy/distributions/collection/Laplace.py
@@ -19,3 +19,12 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.laplace)
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append("loc={loc}")
+        if self.parameters["scale"] != 1.0:
+            s.append("scale={scale}")
+        s = ", ".join(s)
+        return "Laplace(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Levy.py b/src/UQpy/distributions/collection/Levy.py
index b90c4ce10..d55e86951 100644
--- a/src/UQpy/distributions/collection/Levy.py
+++ b/src/UQpy/distributions/collection/Levy.py
@@ -19,3 +19,12 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.levy)
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append("loc={loc}")
+        if self.parameters["scale"] != 1.0:
+            s.append("scale={scale}")
+        s = ", ".join(s)
+        return "Levy(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Logistic.py b/src/UQpy/distributions/collection/Logistic.py
index af3bf69c5..8d4b02cc9 100644
--- a/src/UQpy/distributions/collection/Logistic.py
+++ b/src/UQpy/distributions/collection/Logistic.py
@@ -19,3 +19,12 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.logistic)
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append("loc={loc}")
+        if self.parameters["scale"] != 1.0:
+            s.append("scale={scale}")
+        s = ", ".join(s)
+        return "Logistic(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Lognormal.py b/src/UQpy/distributions/collection/Lognormal.py
index 98b72746a..11cc82b54 100644
--- a/src/UQpy/distributions/collection/Lognormal.py
+++ b/src/UQpy/distributions/collection/Lognormal.py
@@ -25,3 +25,11 @@ def __init__(
             s=s, loc=loc, scale=scale, ordered_parameters=("s", "loc", "scale")
         )
         self._construct_from_scipy(scipy_name=stats.lognorm)
+
+    def __repr__(self):
+        s = "{s}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "Lognormal(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Maxwell.py b/src/UQpy/distributions/collection/Maxwell.py
index 444ea4ebd..5b179600b 100644
--- a/src/UQpy/distributions/collection/Maxwell.py
+++ b/src/UQpy/distributions/collection/Maxwell.py
@@ -19,3 +19,12 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.maxwell)
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append("loc={loc}")
+        if self.parameters["scale"] != 1.0:
+            s.append("scale={scale}")
+        s = ", ".join(s)
+        return "Maxwell(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Multinomial.py b/src/UQpy/distributions/collection/Multinomial.py
index 77cd43a50..6f4ebc309 100644
--- a/src/UQpy/distributions/collection/Multinomial.py
+++ b/src/UQpy/distributions/collection/Multinomial.py
@@ -58,3 +58,6 @@ def moments(self, moments2return="mv"):
             return mean, cov
         else:
             raise ValueError('UQpy: moments2return must be "m", "v" or "mv".')
+
+    def __repr__(self):
+        return "Multinomial({n}, {p})".format(**self.parameters)
diff --git a/src/UQpy/distributions/collection/MultivariateNormal.py b/src/UQpy/distributions/collection/MultivariateNormal.py
index f6fe1d4b1..7d8098250 100644
--- a/src/UQpy/distributions/collection/MultivariateNormal.py
+++ b/src/UQpy/distributions/collection/MultivariateNormal.py
@@ -25,8 +25,13 @@ def __init__(
             if isinstance(cov, (int, float)):
                 pass
             else:
-                if not (len(np.array(cov).shape) in [1, 2] and all(sh == len(mean) for sh in np.array(cov).shape)):
-                    raise ValueError("Input covariance must be a float or ndarray of appropriate dimensions.")
+                if not (
+                    len(np.array(cov).shape) in [1, 2]
+                    and all(sh == len(mean) for sh in np.array(cov).shape)
+                ):
+                    raise ValueError(
+                        "Input covariance must be a float or ndarray of appropriate dimensions."
+                    )
         super().__init__(mean=mean, cov=cov, ordered_parameters=["mean", "cov"])
 
     def cdf(self, x):
@@ -44,8 +49,9 @@ def log_pdf(self, x):
     def rvs(self, nsamples=1, random_state=None):
         if not (isinstance(nsamples, int) and nsamples >= 1):
             raise ValueError("Input nsamples must be an integer > 0.")
-        return stats.multivariate_normal.rvs(size=nsamples, random_state=random_state, **self.parameters
-                                             ).reshape((nsamples, -1))
+        return stats.multivariate_normal.rvs(
+            size=nsamples, random_state=random_state, **self.parameters
+        ).reshape((nsamples, -1))
 
     def fit(self, data):
         data = self.check_x_dimension(data)
@@ -65,3 +71,12 @@ def moments(self, moments2return="mv"):
             return self.get_parameters()["mean"], self.get_parameters()["cov"]
         else:
             raise ValueError('UQpy: moments2return must be "m", "v" or "mv".')
+
+    def __repr__(self):
+        s = "{mean}"
+        is_float_or_int = isinstance(self.parameters["cov"], (int, float))
+        if is_float_or_int:
+            if self.parameters["cov"] == 1.0:
+                return "MultivariateNormal(" + s.format(**self.parameters) + ")"
+        s += ", cov={cov}"
+        return "MultivariateNormal(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Normal.py b/src/UQpy/distributions/collection/Normal.py
index 5c4403909..553961285 100644
--- a/src/UQpy/distributions/collection/Normal.py
+++ b/src/UQpy/distributions/collection/Normal.py
@@ -1,12 +1,14 @@
 from typing import Union
+import numpy as np
 import scipy.stats as stats
 from beartype import beartype
+from scipy.special import erf, erfinv
 from UQpy.distributions.baseclass import DistributionContinuous1D
+from UQpy.utilities.ValidationTypes import NumericArrayLike
 
 
+@beartype
 class Normal(DistributionContinuous1D):
-
-    @beartype
     def __init__(
         self, loc: Union[None, float, int] = 0.0, scale: Union[None, float, int] = 1.0
     ):
@@ -17,3 +19,61 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.norm)
+        self.pdf = self.__probability_density_function
+        self.cdf = self.__cumulative_distribution_function
+        self.icdf = self.__inverse_cumulative_distribution_function
+
+    def __probability_density_function(
+        self, x: NumericArrayLike
+    ) -> Union[float, np.ndarray]:
+        """Probability density function for normal distribution
+
+        :param x: Points to evaluate pdf at
+        :return: PDF of ``x`` as defined by :math:`f_x(x)`
+        """
+        x_array = np.atleast_1d(x)
+        mean = self.parameters["loc"]
+        standard_deviation = self.parameters["scale"]
+        normalizing_constant = 1 / (standard_deviation * np.sqrt(2 * np.pi))
+        pdf = normalizing_constant * np.exp(
+            -0.5 * ((x_array - mean) / standard_deviation) ** 2
+        )
+        return pdf
+
+    def __cumulative_distribution_function(
+        self, x: NumericArrayLike
+    ) -> Union[float, np.ndarray]:
+        """Cumulative distribution function for the normal distribution defined with the error function
+
+        :param x: Points to evaluate cdf at
+        :return: CDF of ``x`` as defined by :math:`F_X(x)`
+        """
+        x_array = np.atleast_1d(x)
+        mean = self.parameters["loc"]
+        standard_deviation = self.parameters["scale"]
+        erf_input = (x_array - mean) / (standard_deviation * np.sqrt(2.0))
+        cdf = (1.0 + erf(erf_input)) / 2.0
+        return cdf
+
+    def __inverse_cumulative_distribution_function(
+        self, y: NumericArrayLike
+    ) -> Union[float, np.ndarray]:
+        """Compute the inverse CDF for the normal distribution with the inverse error function
+
+        :param y:
+        :return:
+        """
+        y_array = np.atleast_1d(y)
+        mean = self.parameters["loc"]
+        standard_deviation = self.parameters["scale"]
+        normalized_icdf = erfinv((2 * y_array) - 1)
+        icdf = (normalized_icdf * standard_deviation * np.sqrt(2.0)) + mean
+        return icdf
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append(f"loc={self.parameters['loc']}")
+        if self.parameters["scale"] != 1.0:
+            s.append(f"scale={self.parameters['scale']}")
+        return "Normal(" + ", ".join(s) + ")"
diff --git a/src/UQpy/distributions/collection/Pareto.py b/src/UQpy/distributions/collection/Pareto.py
index efd6a4f44..28f9494e8 100644
--- a/src/UQpy/distributions/collection/Pareto.py
+++ b/src/UQpy/distributions/collection/Pareto.py
@@ -24,3 +24,11 @@ def __init__(
             b=b, loc=loc, scale=scale, ordered_parameters=("b", "loc", "scale")
         )
         self._construct_from_scipy(scipy_name=stats.pareto)
+
+    def __repr__(self):
+        s = "{b}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "Pareto(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Poisson.py b/src/UQpy/distributions/collection/Poisson.py
index 68027cfb8..b0f4163b9 100644
--- a/src/UQpy/distributions/collection/Poisson.py
+++ b/src/UQpy/distributions/collection/Poisson.py
@@ -17,3 +17,9 @@ def __init__(self, mu: Union[None, float, int], loc: Union[None, float, int] = 0
         """
         super().__init__(mu=mu, loc=loc, ordered_parameters=("mu", "loc"))
         self._construct_from_scipy(scipy_name=stats.poisson)
+
+    def __repr__(self):
+        s = "{mu}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        return "Poisson(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Rayleigh.py b/src/UQpy/distributions/collection/Rayleigh.py
index c66b7d255..42ea3453c 100644
--- a/src/UQpy/distributions/collection/Rayleigh.py
+++ b/src/UQpy/distributions/collection/Rayleigh.py
@@ -18,3 +18,12 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.rayleigh)
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append("loc={loc}")
+        if self.parameters["scale"] != 1.0:
+            s.append("scale={scale}")
+        s = ", ".join(s)
+        return "Rayleigh(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Triangular.py b/src/UQpy/distributions/collection/Triangular.py
new file mode 100644
index 000000000..c7b007867
--- /dev/null
+++ b/src/UQpy/distributions/collection/Triangular.py
@@ -0,0 +1,26 @@
+from scipy import stats
+from UQpy.distributions.baseclass import DistributionContinuous1D
+
+
+class Triangular(DistributionContinuous1D):
+
+    def __init__(self, c: float, loc: float = 0.0, scale: float = 1.0):
+        """
+
+        :param c: Shape parameter between :math:`0 \leq c \leq 1`
+        :param loc: The non-zero part distribution starts at ``loc``. Default: 0.0
+        :param scale: The width of the non-zero part of the distribution. Default: 1.0
+        """
+        super().__init__(
+            c=c, loc=loc, scale=scale, ordered_parameters=("c", "loc", "scale")
+        )
+        self._construct_from_scipy(scipy_name=stats.triang)
+
+    def __repr__(self):
+        s = "c={c}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        s = s.format(**self.parameters)
+        return f"Triangular({s})"
diff --git a/src/UQpy/distributions/collection/TruncatedNormal.py b/src/UQpy/distributions/collection/TruncatedNormal.py
index c862ec469..724d34853 100644
--- a/src/UQpy/distributions/collection/TruncatedNormal.py
+++ b/src/UQpy/distributions/collection/TruncatedNormal.py
@@ -31,3 +31,11 @@ def __init__(
             ordered_parameters=("a", "b", "loc", "scale"),
         )
         self._construct_from_scipy(scipy_name=stats.truncnorm)
+
+    def __repr__(self):
+        s = "{a}, {b}"
+        if self.parameters["loc"] != 0.0:
+            s += ", loc={loc}"
+        if self.parameters["scale"] != 1.0:
+            s += ", scale={scale}"
+        return "TruncatedNormal(" + s.format(**self.parameters) + ")"
diff --git a/src/UQpy/distributions/collection/Uniform.py b/src/UQpy/distributions/collection/Uniform.py
index bc630f218..fcfbd816e 100644
--- a/src/UQpy/distributions/collection/Uniform.py
+++ b/src/UQpy/distributions/collection/Uniform.py
@@ -1,13 +1,13 @@
-from typing import Union
-
+import numpy as np
 import scipy.stats as stats
 from beartype import beartype
-
+from typing import Union
 from UQpy.distributions.baseclass import DistributionContinuous1D
+from UQpy.utilities.ValidationTypes import NumericArrayLike
 
 
+@beartype
 class Uniform(DistributionContinuous1D):
-    @beartype
     def __init__(
         self, loc: Union[None, float, int] = 0.0, scale: Union[None, float, int] = 1.0
     ):
@@ -18,3 +18,69 @@ def __init__(
         """
         super().__init__(loc=loc, scale=scale, ordered_parameters=("loc", "scale"))
         self._construct_from_scipy(scipy_name=stats.uniform)
+
+        self.pdf = self.__probability_density_function
+        self.cdf = self.__cumulative_distribution_function
+        self.icdf = self.__inverse_cumulative_distribution_function
+
+    def __probability_density_function(
+        self, x: NumericArrayLike
+    ) -> Union[int, float, np.ndarray]:
+        """Probability Density Function for the uniform distribution
+
+        :param x: Points at which to evaluate the probability density function
+        :return: pdf at all points in x
+        """
+        x_array = np.atleast_1d(x)
+        loc = self.parameters["loc"]
+        scale = self.parameters["scale"]
+        mask_zero_density = (x_array < loc) | (loc + scale < x_array)
+        mask = (loc <= x_array) & (x_array <= loc + scale)
+        pdf = np.full_like(x_array, np.nan)
+        pdf[mask_zero_density] = 0.0
+        pdf[mask] = 1 / scale
+        return pdf
+
+    def __cumulative_distribution_function(
+        self, x: NumericArrayLike
+    ) -> Union[int, float, np.ndarray]:
+        """Cumulative Distribution Function for the Uniform Distribution
+
+        :param x: Points at which to evaluate the cumulative distribution function
+        :return:  cdf of ``x`` as defined by :math:`F_X(x)`
+        """
+        x_array = np.atleast_1d(x)
+        loc = self.parameters["loc"]
+        scale = self.parameters["scale"]
+        cdf = np.full_like(x_array, np.nan)
+        lower_mask = x_array <= loc
+        middle_mask = (loc < x_array) & (x_array < loc + scale)
+        upper_mask = loc + scale <= x_array
+        cdf[lower_mask] = 0.0
+        cdf[middle_mask] = (x_array[middle_mask] - loc) / scale
+        cdf[upper_mask] = 1.0
+        return cdf
+
+    def __inverse_cumulative_distribution_function(
+        self, x: NumericArrayLike
+    ) -> Union[int, float, np.ndarray]:
+        """Inverse cumulative distribution function for uniform distribution
+
+        :param x: Point at which to evaluate the inverse cumulative distribution function
+        :return: inverse cdf of ``x`` as defined by :math:`F^{-1}_X(x)`
+        """
+        x_array = np.atleast_1d(x)
+        loc = self.parameters["loc"]
+        scale = self.parameters["scale"]
+        icdf = np.full_like(x_array, np.nan)
+        mask = (0 <= x_array) & (x_array <= 1)
+        icdf[mask] = loc + (x_array[mask] * scale)
+        return icdf
+
+    def __repr__(self):
+        s = []
+        if self.parameters["loc"] != 0.0:
+            s.append(f"loc={self.parameters['loc']}")
+        if self.parameters["scale"] != 1.0:
+            s.append(f"scale={self.parameters['scale']}")
+        return "Uniform(" + ", ".join(s) + ")"
diff --git a/src/UQpy/distributions/collection/__init__.py b/src/UQpy/distributions/collection/__init__.py
index 1c4260852..647b7b4e5 100644
--- a/src/UQpy/distributions/collection/__init__.py
+++ b/src/UQpy/distributions/collection/__init__.py
@@ -18,6 +18,7 @@
 from UQpy.distributions.collection.Pareto import Pareto
 from UQpy.distributions.collection.Poisson import Poisson
 from UQpy.distributions.collection.Rayleigh import Rayleigh
+from UQpy.distributions.collection.Triangular import Triangular
 from UQpy.distributions.collection.TruncatedNormal import TruncatedNormal
 from UQpy.distributions.collection.Uniform import Uniform
 from UQpy.distributions.collection.JointIndependent import JointIndependent
diff --git a/src/UQpy/distributions/copulas/Clayton.py b/src/UQpy/distributions/copulas/Clayton.py
index 4cee7922a..640f56238 100644
--- a/src/UQpy/distributions/copulas/Clayton.py
+++ b/src/UQpy/distributions/copulas/Clayton.py
@@ -41,3 +41,6 @@ def extract_data(self, unit_uniform_samples: Numpy2DFloatArray):
         v = unit_uniform_samples[:, 1]
         theta = self.parameters["theta"]
         return theta, u, v
+
+    def __repr__(self):
+        return f"Clayton({self.parameters['theta']})"
diff --git a/src/UQpy/distributions/copulas/Frank.py b/src/UQpy/distributions/copulas/Frank.py
index f633d1e7c..f617cd306 100644
--- a/src/UQpy/distributions/copulas/Frank.py
+++ b/src/UQpy/distributions/copulas/Frank.py
@@ -32,7 +32,11 @@ def evaluate_cdf(self, unit_uniform_samples: Numpy2DFloatArray) -> numpy.ndarray
         :return: Values of the cdf.
         """
         theta, u, v = self.extract_data(unit_uniform_samples)
-        tmp_ratio = ((np.exp(-theta * u) - 1.0) * (np.exp(-theta * v) - 1.0) / (np.exp(-theta) - 1.0))
+        tmp_ratio = (
+            (np.exp(-theta * u) - 1.0)
+            * (np.exp(-theta * v) - 1.0)
+            / (np.exp(-theta) - 1.0)
+        )
         cdf_val = -1.0 / theta * np.log(1.0 + tmp_ratio)
         return cdf_val
 
@@ -41,3 +45,6 @@ def extract_data(self, unit_uniform_samples: Numpy2DFloatArray):
         v = unit_uniform_samples[:, 1]
         theta = self.parameters["theta"]
         return theta, u, v
+
+    def __repr__(self):
+        return f"Frank({self.parameters['theta']})"
diff --git a/src/UQpy/distributions/copulas/Gumbel.py b/src/UQpy/distributions/copulas/Gumbel.py
index d84c1634e..3174714f2 100644
--- a/src/UQpy/distributions/copulas/Gumbel.py
+++ b/src/UQpy/distributions/copulas/Gumbel.py
@@ -73,3 +73,6 @@ def extract_data(self, unit_uniform_samples):
         v = unit_uniform_samples[:, 1]
         theta = self.parameters["theta"]
         return theta, u, v
+
+    def __repr__(self):
+        return f"Gumbel({self.parameters['theta']})"
diff --git a/src/UQpy/utilities/ValidationTypes.py b/src/UQpy/utilities/ValidationTypes.py
index fbd87a966..9c6ed22c4 100644
--- a/src/UQpy/utilities/ValidationTypes.py
+++ b/src/UQpy/utilities/ValidationTypes.py
@@ -4,10 +4,12 @@
 from typing import Union, Annotated
 
 
+
 RandomStateType = Union[None, int, np.random.RandomState]
 PositiveInteger = Annotated[int, Is[lambda number: number > 0]]
 NonNegativeInteger = Annotated[int, Is[lambda number: number >= 0]]
 PositiveFloat = Annotated[float, Is[lambda number: number > 0]]
+NumericArrayLike = Union[int, float, list, tuple, np.ndarray]
 Numpy2DFloatArray = Annotated[
     np.ndarray,
     Is[lambda array: array.ndim == 2 and np.issubdtype(array.dtype, float)],
diff --git a/tests/unit_tests/distributions/test_copula_repr.py b/tests/unit_tests/distributions/test_copula_repr.py
new file mode 100644
index 000000000..3f6e436bf
--- /dev/null
+++ b/tests/unit_tests/distributions/test_copula_repr.py
@@ -0,0 +1,17 @@
+import numpy as np
+from UQpy.distributions import copulas
+
+
+def test_clayton_repr():
+    clayton = copulas.Clayton(1.0)
+    assert clayton.__repr__() == "Clayton(1.0)"
+
+
+def test_frank_repr():
+    frank = copulas.Frank(2.0)
+    assert frank.__repr__() == "Frank(2.0)"
+
+
+def test_gumbel_repr():
+    gumbel = copulas.Gumbel(np.inf)
+    assert gumbel.__repr__() == "Gumbel(inf)"
diff --git a/tests/unit_tests/distributions/test_distribution_methods.py b/tests/unit_tests/distributions/test_distribution_methods.py
index f1afa352a..34dee4d24 100644
--- a/tests/unit_tests/distributions/test_distribution_methods.py
+++ b/tests/unit_tests/distributions/test_distribution_methods.py
@@ -1,18 +1,19 @@
-from UQpy.distributions import *
 import numpy as np
+from UQpy.distributions import *
+
 
 # Test all functions for one type of continuous distribution: uniform
-dist_continuous = Uniform(loc=1., scale=2.)
+dist_continuous = Uniform(loc=1.0, scale=2.0)
 
 
 def test_get_params():
-    assert dist_continuous.get_parameters()['loc'] == 1.
+    assert dist_continuous.get_parameters()["loc"] == 1.0
 
 
 def test_update_params():
-    dist = Uniform(loc=1., scale=2.)
-    dist.update_parameters(loc=2.)
-    assert dist.get_parameters()['loc'] == 2.
+    dist = Uniform(loc=1.0, scale=2.0)
+    dist.update_parameters(loc=2.0)
+    assert dist.get_parameters()["loc"] == 2.0
 
 
 def test_continuous_pdf():
@@ -38,11 +39,11 @@ def test_continuous_rvs():
 
 def test_continuous_fit():
     dict_fit = Uniform(loc=None, scale=None).fit(data=[1.5, 2.5, 3.5])
-    assert dict_fit == {'loc': 1.5, 'scale': 2.0}
+    assert dict_fit == {"loc": 1.5, "scale": 2.0}
 
 
 def test_continuous_moments():
-    assert dist_continuous.moments(moments2return='m') == 2.
+    assert dist_continuous.moments(moments2return="m") == 2.0
 
 
 # Test all functions for one type of discrete distribution: binomial
@@ -50,34 +51,32 @@ def test_continuous_moments():
 
 
 def test_discrete_pmf():
-    assert np.round(dist_discrete.pmf(x=2.), 3) == 0.205
+    assert np.round(dist_discrete.pmf(x=2.0), 3) == 0.205
 
 
 def test_discrete_cdf():
-    assert np.round(dist_discrete.cdf(x=2.), 3) == 0.942
+    assert np.round(dist_discrete.cdf(x=2.0), 3) == 0.942
 
 
 def test_discrete_log_pmf():
-    assert np.round(dist_discrete.log_pmf(x=2.), 3) == -1.586
+    assert np.round(dist_discrete.log_pmf(x=2.0), 3) == -1.586
 
 
 def test_discrete_icdf():
-    assert dist_discrete.icdf(0.9) == 2.
+    assert dist_discrete.icdf(0.9) == 2.0
 
 
 def test_discrete_rvs():
     samples = dist_discrete.rvs(nsamples=2, random_state=123)
-    assert np.all(np.round(samples, 3) == np.array([1., 0.]).reshape((2, 1)))
+    assert np.all(np.round(samples, 3) == np.array([1.0, 0.0]).reshape((2, 1)))
 
 
 def test_discrete_moments():
-    assert dist_discrete.moments(moments2return='m') == 1.
+    assert dist_discrete.moments(moments2return="m") == 1.0
 
 
 # Test functions for Copula
-
-
 def test_update_params_copula():
-    copula = Gumbel(theta=2.)
-    copula.update_parameters(theta=1.)
-    assert copula.get_parameters()['theta'] == 1.
+    copula = Gumbel(theta=2.0)
+    copula.update_parameters(theta=1.0)
+    assert copula.get_parameters()["theta"] == 1.0
diff --git a/tests/unit_tests/distributions/test_distribution_repr.py b/tests/unit_tests/distributions/test_distribution_repr.py
new file mode 100644
index 000000000..1fcf20207
--- /dev/null
+++ b/tests/unit_tests/distributions/test_distribution_repr.py
@@ -0,0 +1,129 @@
+from UQpy import distributions
+
+
+def test_beta_repr():
+    beta = distributions.Beta(1, 2, 3, 4)
+    assert beta.__repr__() == "Beta(1, 2, loc=3, scale=4)"
+
+
+def test_binomial_repr():
+    binomial = distributions.Binomial(1, 2, 3)
+    assert binomial.__repr__() == "Binomial(1, 2, loc=3)"
+
+
+def test_cauchy_repr():
+    cauchy = distributions.Cauchy(1, 2)
+    assert cauchy.__repr__() == "Cauchy(loc=1, scale=2)"
+
+
+def test_chi_square_repr():
+    chi_square = distributions.ChiSquare(1, 2, 3)
+    assert chi_square.__repr__() == "ChiSquare(1, loc=2, scale=3)"
+
+
+def test_exponential_repr():
+    exponential = distributions.Exponential(1, 2)
+    assert exponential.__repr__() == "Exponential(loc=1, scale=2)"
+
+
+def test_gamma_repr():
+    gamma = distributions.Gamma(1, 2, 3)
+    assert gamma.__repr__() == "Gamma(1, loc=2, scale=3)"
+
+
+def test_generalized_extreme_repr():
+    generalized_extreme = distributions.GeneralizedExtreme(1, 2, 3)
+    assert generalized_extreme.__repr__() == "GeneralizedExtreme(1, loc=2, scale=3)"
+
+
+def test_inverse_gaussian_repr():
+    inverse_gauss = distributions.InverseGauss(1, 2, 3)
+    assert inverse_gauss.__repr__() == "InverseGauss(1, loc=2, scale=3)"
+
+
+def test_joint_copula_repr():
+    joint_copula = distributions.JointCopula(
+        [distributions.Normal(loc=2.0), distributions.Uniform()],
+        copula=distributions.copulas.Frank(1.0),
+    )
+    assert (
+        joint_copula.__repr__()
+        == "JointCopula([Normal(loc=2.0), Uniform()], Frank(1.0))"
+    )
+
+
+def test_joint_independent_repr():
+    joint_independent = distributions.JointIndependent(
+        [distributions.Beta(1, 2), distributions.Cauchy()]
+    )
+    assert joint_independent.__repr__() == "JointIndependent([Beta(1, 2), Cauchy()])"
+
+
+def test_laplace_repr():
+    laplace = distributions.Laplace(1, 2)
+    assert laplace.__repr__() == "Laplace(loc=1, scale=2)"
+
+
+def test_levy_repr():
+    levy = distributions.Levy(1, 2)
+    assert levy.__repr__() == "Levy(loc=1, scale=2)"
+
+
+def test_logistic_repr():
+    logistic = distributions.Logistic(1, 2)
+    assert logistic.__repr__() == "Logistic(loc=1, scale=2)"
+
+
+def test_lognormal_repr():
+    lognormal = distributions.Lognormal(1, 2, 3)
+    assert lognormal.__repr__() == "Lognormal(1, loc=2, scale=3)"
+
+
+def test_maxwell_repr():
+    maxwell = distributions.Maxwell(1, 2)
+    assert maxwell.__repr__() == "Maxwell(loc=1, scale=2)"
+
+
+def test_multinomial_repr():
+    multinomial = distributions.Multinomial(1, [2.0, 3.0])
+    assert multinomial.__repr__() == "Multinomial(1, [2.0, 3.0])"
+
+
+def test_multivariate_normal_repr():
+    multivariate_normal = distributions.MultivariateNormal(
+        [1, 2], [[1.0, 0.4], [0.4, 1.0]]
+    )
+    assert (
+        multivariate_normal.__repr__()
+        == "MultivariateNormal([1, 2], cov=[[1.0, 0.4], [0.4, 1.0]])"
+    )
+
+
+def test_normal_repr():
+    normal = distributions.Normal(1, 2)
+    assert normal.__repr__() == "Normal(loc=1, scale=2)"
+
+
+def test_pareto_repr():
+    pareto = distributions.Pareto(1, 2, 3)
+    assert pareto.__repr__() == "Pareto(1, loc=2, scale=3)"
+
+
+def test_poisson_repr():
+    poisson = distributions.Poisson(1, 2)
+    assert poisson.__repr__() == "Poisson(1, loc=2)"
+
+
+def test_rayleigh_repr():
+    rayleigh = distributions.Rayleigh(1, 2)
+    assert rayleigh.__repr__() == "Rayleigh(loc=1, scale=2)"
+
+
+def test_truncated_normal_repr():
+    truncated_normal = distributions.TruncatedNormal(1, 2, 3, 4)
+    assert truncated_normal.__repr__() == "TruncatedNormal(1, 2, loc=3, scale=4)"
+
+
+def test_uniform_repr():
+    uniform = distributions.Uniform(1, 2)
+    assert uniform.__repr__() == "Uniform(loc=1, scale=2)"
diff --git a/tests/unit_tests/distributions/test_normal.py b/tests/unit_tests/distributions/test_normal.py
new file mode 100644
index 000000000..817fe2ea5
--- /dev/null
+++ b/tests/unit_tests/distributions/test_normal.py
@@ -0,0 +1,81 @@
+import pytest
+import numpy as np
+from scipy import stats
+from UQpy.distributions import Normal
+from hypothesis import given
+from hypothesis.extra.numpy import array_shapes
+
+
+class TestNormal:
+
+    normal = Normal()
+    scipy_normal = stats.norm()
+
+    def test_pdf_nan(self):
+        """Consistent with scipy, pdf(NaN)=NaN"""
+        assert np.isnan(self.normal.pdf(np.nan))
+
+    def test_pdf_infinity(self):
+        """Consistent with scipy pdf(inf)=0.0, pdf(-inf)=0.0"""
+        assert np.isclose(self.normal.pdf(np.inf), 0.0)
+        assert np.isclose(self.normal.pdf(-np.inf), 0.0)
+
+    @given(array_shapes(min_dims=1, min_side=1))
+    def test_pdf_shape(self, shape):
+        """Test the output array matches the shape of the input array"""
+        x = np.zeros(shape)
+        assert x.shape == self.normal.pdf(x).shape
+
+    def test_pdf_values(self):
+        """Test if UQpy pdf matches SciPy pdf for x in [-10, 10]"""
+        x = np.linspace(-10, 10, num=1_000)
+        assert np.allclose(self.normal.pdf(x), self.scipy_normal.pdf(x))
+
+    def test_cdf_nan(self):
+        """Consistent with scipy, cdf(NaN)=NaN"""
+        assert np.isnan(self.normal.cdf(np.nan))
+
+    @pytest.mark.parametrize("test_input,expected", [(np.inf, 1.0), (-np.inf, 0.0)])
+    def test_cdf_infinity(self, test_input, expected):
+        """Consistent with scipy cdf(inf)=1.0, cdf(-inf)=0.0"""
+        assert self.normal.cdf(test_input) == expected
+
+    @given(array_shapes(min_dims=1, min_side=1))
+    def test_cdf_shape(self, shape):
+        """Test if output array matches the shape of the input array"""
+        x = np.zeros(shape)
+        assert x.shape == self.normal.cdf(x).shape
+
+    def test_cdf_values(self):
+        """Test if UQpy cdf matches scipy cdf for x in [-10, 10]"""
+        x = np.linspace(-10, 10, num=100)
+        assert np.allclose(self.normal.cdf(x), self.scipy_normal.cdf(x))
+
+    def test_icdf_nan(self):
+        """Consistent with scipy, icdf(NaN) = NaN"""
+        assert np.isnan(self.normal.icdf(np.nan))
+
+    def test_icdf_infinity(self):
+        """Consistent with scipy icdf(inf)=NaN and icdf(-inf)=NaN"""
+        assert np.isnan(self.normal.icdf(np.inf))
+        assert np.isnan(self.normal.icdf(-np.inf))
+
+    def test_icdf_zero_one(self):
+        """Consistent with scipy, icdf(0)=-inf and icdf(1)=inf"""
+        assert self.normal.icdf(0) == -np.inf
+        assert self.normal.icdf(1) == np.inf
+
+    def test_icdf_values(self):
+        """Test if UQpy icdf matches scipy ppf for y in [-0.1, 1.1]. Note outside [0, 1] both functions return NaN"""
+        y = np.linspace(-0.1, 1.1, num=100)
+        assert np.allclose(self.normal.icdf(y), self.scipy_normal.ppf(y), equal_nan=True)
+
+    def test_cdf_icdf_reconstruction(self):
+        """Reconstruct x as x = icdf(cdf(x))
+        Note: For +/- 7 standard deviations, UQpy and SciPy accurately reconstruct x.
+        At 8 standard deviations, both UQpy and scipy.stats.norm() begin to divergence from the correct answer
+        """
+        x = np.linspace(-7, 7, num=100)
+        y = self.normal.cdf(x)
+        x_reconstruction = self.normal.icdf(y)
+        assert np.allclose(x, x_reconstruction)
diff --git a/tests/unit_tests/distributions/test_uniform.py b/tests/unit_tests/distributions/test_uniform.py
new file mode 100644
index 000000000..0eea47206
--- /dev/null
+++ b/tests/unit_tests/distributions/test_uniform.py
@@ -0,0 +1,81 @@
+import pytest
+import numpy as np
+from scipy import stats
+from UQpy.distributions import Uniform
+from hypothesis import given
+from hypothesis.extra.numpy import array_shapes
+
+
+class TestUniform:
+
+    uniform = Uniform()
+    scipy_uniform = stats.uniform()
+
+    def test_pdf_nan(self):
+        """Consistent with scipy, pdf(NaN)=NaN"""
+        assert np.isnan(self.uniform.pdf(np.nan))
+
+    def test_pdf_infinity(self):
+        """Consistent with scipy, pdf(inf)=0 and pdf(-inf)=0"""
+        assert np.isclose(self.uniform.pdf(np.inf), 0.0)
+        assert np.isclose(self.uniform.pdf(-np.inf), 0.0)
+
+    @given(array_shapes(min_dims=1, min_side=1))
+    def test_pdf_shape(self, shape):
+        """Test if the output array matches the shape of the input array"""
+        x = np.zeros(shape)
+        assert x.shape == self.uniform.pdf(x).shape
+
+    def test_pdf_values(self):
+        """Test if UQpy pdf matches scipy pdf for x in [-1, 2]"""
+        x = np.linspace(-1, 2, num=100)
+        assert np.allclose(self.uniform.pdf(x), self.scipy_uniform.pdf(x))
+
+    def test_cdf_nan(self):
+        """Consistent with scipy cdf(NaN)=NaN"""
+        assert np.isnan(self.uniform.cdf(np.nan))
+
+    @pytest.mark.parametrize("test_input,expected", [(np.inf, 1.0), (-np.inf, 0.0)])
+    def test_cdf_infinity(self, test_input, expected):
+        """Consistent with scipy cdf(inf)=1.0, cdf(-inf)=0.0"""
+        assert self.uniform.cdf(test_input) == expected
+
+    @given(array_shapes(min_dims=1, min_side=1))
+    def test_cdf_shape(self, shape):
+        """Test if the output array matches the shape of the input array"""
+        x = np.zeros(shape)
+        assert x.shape == self.uniform.cdf(x).shape
+
+    def test_cdf_values(self):
+        """Test if UQpy cdf matches scipy cdf on x in [-1, 2]"""
+        x = np.linspace(-1, 2, num=100)
+        assert np.allclose(self.uniform.cdf(x), self.scipy_uniform.cdf(x))
+
+    def test_icdf_nan(self):
+        """Consistent with scipy icdf(NaN)=NaN"""
+        assert np.isnan(self.uniform.icdf(np.nan))
+
+    def test_icdf_infinity(self):
+        """Consistent with scipy icdf(inf)=NaN and icdf(-inf)=NaN"""
+        assert np.isnan(self.uniform.icdf(np.inf))
+        assert np.isnan(self.uniform.icdf(-np.inf))
+
+    @given(array_shapes(min_dims=1, min_side=1))
+    def test_icdf_shape(self, shape):
+        """Test if the output array has the same shape as the input array"""
+        x = np.zeros(shape)
+        assert x.shape == self.uniform.icdf(x).shape
+
+    def test_icdf_values(self):
+        """Test if UQpy icdf matches scipy ppf for y in [-0.1, 1.1]. Note outside of [0, 1] both return NaN"""
+        y = np.linspace(-0.1, 1.1, num=100)
+        assert np.allclose(
+            self.uniform.icdf(y), self.scipy_uniform.ppf(y), equal_nan=True
+        )
+
+    def test_cdf_icdf_reconstruction(self):
+        """Reconstruct x as x = icdf(cdf(x)). Note that this only works where cdf(x) is invertible on [0,1]"""
+        x = np.linspace(0, 1, num=100)
+        y = self.uniform.cdf(x)
+        x_reconstruction = self.uniform.icdf(y)
+        assert np.allclose(x, x_reconstruction)