From 7f40006a0e44c6cc8187d87ea3769f34072467c2 Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Thu, 24 Sep 2020 21:39:06 -0400
Subject: [PATCH 01/14] HyperLaplace distribution implemented

---
 copulas/univariate/__init__.py      |   4 +-
 copulas/univariate/hyper_laplace.py | 282 ++++++++++++++++++++++++++++
 2 files changed, 285 insertions(+), 1 deletion(-)
 create mode 100644 copulas/univariate/hyper_laplace.py

diff --git a/copulas/univariate/__init__.py b/copulas/univariate/__init__.py
index 954e5ee1..399bf3d3 100644
--- a/copulas/univariate/__init__.py
+++ b/copulas/univariate/__init__.py
@@ -7,6 +7,7 @@
 from copulas.univariate.student_t import StudentTUnivariate
 from copulas.univariate.truncated_gaussian import TruncatedGaussian
 from copulas.univariate.uniform import UniformUnivariate
+from copulas.univariate.hyper_laplace import HyperLaplace
 
 __all__ = (
     'BetaUnivariate',
@@ -19,5 +20,6 @@
     'ParametricType',
     'BoundedType',
     'UniformUnivariate',
-    'LogLaplace'
+    'LogLaplace',
+    'HyperLaplace'
 )
diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
new file mode 100644
index 00000000..c00f5409
--- /dev/null
+++ b/copulas/univariate/hyper_laplace.py
@@ -0,0 +1,282 @@
+import numpy as np
+from scipy.optimize import fmin_slsqp
+import scipy.stats
+from random import random
+
+from copulas import EPSILON, store_args
+from copulas.univariate.base import BoundedType, ParametricType, Univariate
+from scipy.special import gamma, psi
+
+def gamma_d(x):
+    """
+    Compute the derivative of gamma function
+
+    Arguments:
+            x (float)
+
+        Returns:
+            float:
+                The derivative of gamma function at x
+    """
+    return gamma(x)*psi(x)
+
+def g3g1_over_g22(x):
+    """
+    Commute the value of Gamma(3x)Gamma(x)/Gamma^2(2x)
+
+    Arguments:
+            x (float)
+
+        Returns:
+            float:
+                The value of Gamma(3x)Gamma(x)/Gamma^2(2x)
+    """
+    return gamma(3*x)*gamma(x)/gamma(2*x)**2
+
+def g3g1_over_g22_d(x):
+    """
+    Compute the derivative of Gamma(3x)Gamma(x)/Gamma^2(2x)
+
+    Arguments:
+            x (float)
+
+        Returns:
+            float:
+                The derivative of Gamma(3x)Gamma(x)/Gamma^2(2x) at x
+    """
+    x1 = gamma(x)
+    x2 = gamma(2*x)
+    x3 = gamma(3*x)
+    
+    x3_d = 3*gamma_d(3*x)*x1*(x2**2)
+    x1_d = x3*gamma_d(x)*(x2**2)
+    x2_d = 4*x3*x1*gamma_d(2*x)*x2
+    return (x3_d + x1_d - x2_d)/(x2**4)
+
+def solve_gamma_equation(y, num_iter = 10, upper_threshold = 2, lower_threshold = 1/2):
+    """
+    Solve for the equation Gamma(3x)Gamma(x)/Gamma^2(2x) = y using Newton's method
+
+    Arguments:
+            y (float)
+
+            num_iter (int):
+                number of iterations when using Newton's method
+            
+            upper_threshold (float):
+                Maximum value of x_new/x_old allowed in Newton's method. Should be > 1
+
+            lower_threshold (float):
+                Minimum value of x_new/x_old allowed in Newton's method. Should be in (0,1)
+
+        Returns:
+            float:
+                The solution to Gamma(3x)Gamma(x)/Gamma^2(2x) = y
+    """
+    x = 1
+    for i in range(num_iter):
+        x_new = x + (y - g3g1_over_g22(x))/g3g1_over_g22_d(x)
+        if x_new > upper_threshold*x:
+            x = upper_threshold*x
+        elif x_new < lower_threshold*x:
+            x = lower_threshold*x
+        else:
+            x = x_new
+    return x
+
+
+class HyperLaplace(Univariate):
+    """
+    An HyperLaplace model object, implemented via a wrapper around scipy.stats.gamma
+
+    Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html
+
+    Math derivation: HyperLaplace(k, alpha) = (gamma(loc = 0, scale = 1/k, a = 1/alpha))**(1/alpha) * Unif({-1,1})
+    """
+
+    MODEL_CLASS = scipy.stats.gamma
+
+    _params = None
+    _model = None
+    
+    alpha = None
+    k = None
+
+    def probability_density(self, X):
+        """Compute the probability density for each point in X.
+
+        Arguments:
+            X (numpy.ndarray):
+                Values for which the probability density will be computed.
+                It must have shape (n, 1).
+
+        Returns:
+            numpy.ndarray:
+                Probability density values for points in X.
+
+        Raises:
+            NotFittedError:
+                if the model is not fitted.
+        """
+        self.check_fit()
+        return self._model.pdf(abs(X)**self.alpha)*0.5*alpha*abs(X)**(self.alpha - 1)
+
+    def log_probability_density(self, X):
+        """Compute the log of the probability density for each point in X.
+
+        Arguments:
+            X (numpy.ndarray):
+                Values for which the log probability density will be computed.
+                It must have shape (n, 1).
+
+        Returns:
+            numpy.ndarray:
+                Log probability density values for points in X.
+
+        Raises:
+            NotFittedError:
+                if the model is not fitted.
+        """
+        self.check_fit()
+
+        return np.log(self.probability_density(X))
+
+    def cumulative_distribution(self, X):
+        """Compute the cumulative distribution value for each point in X.
+
+        Arguments:
+            X (numpy.ndarray):
+                Values for which the cumulative distribution will be computed.
+                It must have shape (n, 1).
+
+        Returns:
+            numpy.ndarray:
+                Cumulative distribution values for points in X.
+
+        Raises:
+            NotFittedError:
+                if the model is not fitted.
+        """
+        self.check_fit()
+        gamma_cdf = self._model.cdf(abs(X)**self.alpha)
+        for i in range(len(X)):
+            cd, x = gamma_cdf[i], X[i]
+            if x > 0:
+                gamma_cdf[i] = (cd + 1)/2
+            else:
+                gamma_cdf[i] = (1-cd)/2
+
+        return gamma_cdf
+
+    def percent_point(self, U):
+        """Compute the inverse cumulative distribution value for each point in U.
+
+        Arguments:
+            U (numpy.ndarray):
+                Values for which the cumulative distribution will be computed.
+                It must have shape (n, 1) and values must be in [0,1].
+
+        Returns:
+            numpy.ndarray:
+                Inverse cumulative distribution values for points in U.
+
+        Raises:
+            NotFittedError:
+                if the model is not fitted.
+        """
+        self.check_fit()
+        gamma_ppf = self._model.ppf(abs(U-0.5)*2)**(1/self.alpha)
+        for i in range(len(U)):
+            pp, u = gamma_ppf[i], U[i]
+            if u < 0.5:
+                gamma_ppf[i] = -pp
+
+        return gamma_ppf
+
+    def sample(self, n_samples=1):
+        """Sample values from this model.
+
+        Argument:
+            n_samples (int):
+                Number of values to sample
+
+        Returns:
+            numpy.ndarray:
+                Array of shape (n_samples, 1) with values randomly
+                sampled from this model distribution.
+
+        Raises:
+            NotFittedError:
+                if the model is not fitted.
+        """
+        self.check_fit()
+        gamma_rvs = self._model.rvs(n_samples)**(1/self.alpha)
+        for i in range(n_samples):
+            if random() > 0.5:
+                gamma_rvs[i] = -gamma_rvs[i]
+        return gamma_rvs
+
+    def _fit(self, X):
+        """Fit the model to a non-constant random variable.
+
+        This fitting method is implemented by matching the theoritical mean and variance with the emperical mean and variance.
+
+        Arguments:
+            X (numpy.ndarray):
+                Values of the random variable. It must have shape (n, 1).
+        """
+        mean = np.mean(abs(X))
+        square = np.mean(X**2)
+        a = solve_gamma_equation(square/(mean**2))
+        self.alpha = 1/a
+        self.k = (gamma(2*a)/(gamma(a)*mean))**self.alpha
+
+        self._params = {
+            'loc': 0,
+            'scale': 1/self.k,
+            'a': a
+        }
+
+    def _get_model(self):
+        return self.MODEL_CLASS(**self._params)
+
+    def fit(self, X):
+        """Fit the model to a random variable.
+
+        Arguments:
+            X (numpy.ndarray):
+                Values of the random variable. It must have shape (n, 1).
+        """
+        if self._check_constant_value(X):
+            self._fit_constant(X)
+        else:
+            self._fit(X)
+            self._model = self._get_model()
+
+        self.fitted = True
+
+    def _get_params(self):
+        """Return attributes from self._model to serialize.
+
+        Must be implemented in all the subclasses.
+
+        Returns:
+            dict:
+                Parameters to recreate self._model in its current fit status.
+        """
+        return self._params.copy()
+
+    def _set_params(self, params):
+        """Set the parameters of this univariate.
+
+        Args:
+            params (dict):
+                Parameters to recreate this instance.
+        """
+        self._params = params.copy()
+        self.k = 1/self._params['scale']
+        self.alpha = 1/self._params['a']
+        if self._is_constant():
+            self._replace_constant_methods()
+        else:
+            self._model = self._get_model()
\ No newline at end of file

From 5e8e18f56456f1cfa2fe7d2a670aa51d8d3c612f Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Fri, 25 Sep 2020 13:16:56 -0400
Subject: [PATCH 02/14] Modified tests to include hyper_laplace as a candidate

---
 copulas/univariate/hyper_laplace.py | 14 ++++++++++++++
 tests/unit/univariate/test_base.py  | 10 +++++++---
 2 files changed, 21 insertions(+), 3 deletions(-)

diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index c00f5409..58b31c64 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -94,6 +94,9 @@ class HyperLaplace(Univariate):
     Math derivation: HyperLaplace(k, alpha) = (gamma(loc = 0, scale = 1/k, a = 1/alpha))**(1/alpha) * Unif({-1,1})
     """
 
+    PARAMETRIC = ParametricType.PARAMETRIC
+    BOUNDED = BoundedType.UNBOUNDED
+
     MODEL_CLASS = scipy.stats.gamma
 
     _params = None
@@ -240,6 +243,17 @@ def _fit(self, X):
     def _get_model(self):
         return self.MODEL_CLASS(**self._params)
 
+
+    def _fit_constant(self, X):
+        self._params = {
+            'loc': np.unique(X)[0],
+            'scale': 0,
+            'a': 1
+        }
+
+    def _is_constant(self):
+        return self._params['scale'] == 0
+        
     def fit(self, X):
         """Fit the model to a random variable.
 
diff --git a/tests/unit/univariate/test_base.py b/tests/unit/univariate/test_base.py
index f76c3b0b..16c2f2b7 100644
--- a/tests/unit/univariate/test_base.py
+++ b/tests/unit/univariate/test_base.py
@@ -11,6 +11,7 @@
 from copulas.univariate.student_t import StudentTUnivariate
 from copulas.univariate.truncated_gaussian import TruncatedGaussian
 from copulas.univariate.uniform import UniformUnivariate
+from copulas.univariate.hyper_laplace import HyperLaplace
 from tests import compare_nested_iterables
 
 
@@ -29,7 +30,8 @@ def test__select_candidates(self):
             GammaUnivariate,
             StudentTUnivariate,
             UniformUnivariate,
-            LogLaplace
+            LogLaplace,
+            HyperLaplace
         }
 
     def test__select_candidates_parametric(self):
@@ -44,7 +46,8 @@ def test__select_candidates_parametric(self):
             GammaUnivariate,
             StudentTUnivariate,
             UniformUnivariate,
-            LogLaplace
+            LogLaplace,
+            HyperLaplace
         }
 
     def test__select_candidates_non_parametric(self):
@@ -73,7 +76,8 @@ def test__select_candidates_unbounded(self):
         assert set(candidates) == {
             GaussianKDE,
             GaussianUnivariate,
-            StudentTUnivariate
+            StudentTUnivariate,
+            HyperLaplace
         }
 
     def test__select_candidates_semibounded(self):

From aa394ade4baea85bee25cb0d4ce99c78534e5032 Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Fri, 25 Sep 2020 13:34:19 -0400
Subject: [PATCH 03/14] modified code style

---
 copulas/univariate/hyper_laplace.py | 67 ++++++++++++++---------------
 1 file changed, 32 insertions(+), 35 deletions(-)

diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index 58b31c64..804c93ea 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -1,9 +1,7 @@
 import numpy as np
-from scipy.optimize import fmin_slsqp
 import scipy.stats
 from random import random
 
-from copulas import EPSILON, store_args
 from copulas.univariate.base import BoundedType, ParametricType, Univariate
 from scipy.special import gamma, psi
 
@@ -18,7 +16,7 @@ def gamma_d(x):
             float:
                 The derivative of gamma function at x
     """
-    return gamma(x)*psi(x)
+    return gamma(x) * psi(x)
 
 def g3g1_over_g22(x):
     """
@@ -31,7 +29,7 @@ def g3g1_over_g22(x):
             float:
                 The value of Gamma(3x)Gamma(x)/Gamma^2(2x)
     """
-    return gamma(3*x)*gamma(x)/gamma(2*x)**2
+    return gamma(3 * x) * gamma(x) / gamma(2 * x) ** 2
 
 def g3g1_over_g22_d(x):
     """
@@ -45,15 +43,15 @@ def g3g1_over_g22_d(x):
                 The derivative of Gamma(3x)Gamma(x)/Gamma^2(2x) at x
     """
     x1 = gamma(x)
-    x2 = gamma(2*x)
-    x3 = gamma(3*x)
+    x2 = gamma(2 * x)
+    x3 = gamma(3 * x)
     
-    x3_d = 3*gamma_d(3*x)*x1*(x2**2)
-    x1_d = x3*gamma_d(x)*(x2**2)
-    x2_d = 4*x3*x1*gamma_d(2*x)*x2
-    return (x3_d + x1_d - x2_d)/(x2**4)
+    x3_d = 3 * gamma_d(3 * x) * x1 * (x2 ** 2)
+    x1_d = x3 * gamma_d(x) * (x2 ** 2)
+    x2_d = 4 * x3 * x1 * gamma_d(2 * x) * x2
+    return (x3_d + x1_d - x2_d) / (x2 ** 4)
 
-def solve_gamma_equation(y, num_iter = 10, upper_threshold = 2, lower_threshold = 1/2):
+def solve_gamma_equation(y, num_iter = 10, upper_threshold = 2, lower_threshold = 0.5):
     """
     Solve for the equation Gamma(3x)Gamma(x)/Gamma^2(2x) = y using Newton's method
 
@@ -75,11 +73,11 @@ def solve_gamma_equation(y, num_iter = 10, upper_threshold = 2, lower_threshold
     """
     x = 1
     for i in range(num_iter):
-        x_new = x + (y - g3g1_over_g22(x))/g3g1_over_g22_d(x)
-        if x_new > upper_threshold*x:
-            x = upper_threshold*x
-        elif x_new < lower_threshold*x:
-            x = lower_threshold*x
+        x_new = x + (y - g3g1_over_g22(x)) / g3g1_over_g22_d(x)
+        if x_new > upper_threshold * x:
+            x = upper_threshold * x
+        elif x_new < lower_threshold * x:
+            x = lower_threshold * x
         else:
             x = x_new
     return x
@@ -122,7 +120,7 @@ def probability_density(self, X):
                 if the model is not fitted.
         """
         self.check_fit()
-        return self._model.pdf(abs(X)**self.alpha)*0.5*alpha*abs(X)**(self.alpha - 1)
+        return self._model.pdf(abs(X) ** self.alpha) * 0.5 * alpha * abs(X) ** (self.alpha - 1)
 
     def log_probability_density(self, X):
         """Compute the log of the probability density for each point in X.
@@ -161,13 +159,13 @@ def cumulative_distribution(self, X):
                 if the model is not fitted.
         """
         self.check_fit()
-        gamma_cdf = self._model.cdf(abs(X)**self.alpha)
+        gamma_cdf = self._model.cdf(abs(X) ** self.alpha)
         for i in range(len(X)):
             cd, x = gamma_cdf[i], X[i]
             if x > 0:
-                gamma_cdf[i] = (cd + 1)/2
+                gamma_cdf[i] = (cd + 1) / 2
             else:
-                gamma_cdf[i] = (1-cd)/2
+                gamma_cdf[i] = (1 - cd) / 2
 
         return gamma_cdf
 
@@ -188,15 +186,15 @@ def percent_point(self, U):
                 if the model is not fitted.
         """
         self.check_fit()
-        gamma_ppf = self._model.ppf(abs(U-0.5)*2)**(1/self.alpha)
+        gamma_ppf = self._model.ppf(abs(U - 0.5) * 2) ** (1 / self.alpha)
         for i in range(len(U)):
             pp, u = gamma_ppf[i], U[i]
             if u < 0.5:
-                gamma_ppf[i] = -pp
+                gamma_ppf[i] = - pp
 
         return gamma_ppf
 
-    def sample(self, n_samples=1):
+    def sample(self, n_samples = 1):
         """Sample values from this model.
 
         Argument:
@@ -213,10 +211,10 @@ def sample(self, n_samples=1):
                 if the model is not fitted.
         """
         self.check_fit()
-        gamma_rvs = self._model.rvs(n_samples)**(1/self.alpha)
+        gamma_rvs = self._model.rvs(n_samples) ** (1 / self.alpha)
         for i in range(n_samples):
             if random() > 0.5:
-                gamma_rvs[i] = -gamma_rvs[i]
+                gamma_rvs[i] = - gamma_rvs[i]
         return gamma_rvs
 
     def _fit(self, X):
@@ -229,21 +227,20 @@ def _fit(self, X):
                 Values of the random variable. It must have shape (n, 1).
         """
         mean = np.mean(abs(X))
-        square = np.mean(X**2)
-        a = solve_gamma_equation(square/(mean**2))
-        self.alpha = 1/a
-        self.k = (gamma(2*a)/(gamma(a)*mean))**self.alpha
+        square = np.mean(X ** 2)
+        a = solve_gamma_equation(square / (mean **2 ))
+        self.alpha = 1 / a
+        self.k = (gamma(2 * a) / (gamma(a) * mean)) ** self.alpha
 
         self._params = {
             'loc': 0,
-            'scale': 1/self.k,
+            'scale': 1 / self.k,
             'a': a
         }
 
     def _get_model(self):
         return self.MODEL_CLASS(**self._params)
 
-
     def _fit_constant(self, X):
         self._params = {
             'loc': np.unique(X)[0],
@@ -253,7 +250,7 @@ def _fit_constant(self, X):
 
     def _is_constant(self):
         return self._params['scale'] == 0
-        
+
     def fit(self, X):
         """Fit the model to a random variable.
 
@@ -288,9 +285,9 @@ def _set_params(self, params):
                 Parameters to recreate this instance.
         """
         self._params = params.copy()
-        self.k = 1/self._params['scale']
-        self.alpha = 1/self._params['a']
+        self.k = 1 / self._params['scale']
+        self.alpha = 1 / self._params['a']
         if self._is_constant():
             self._replace_constant_methods()
         else:
-            self._model = self._get_model()
\ No newline at end of file
+            self._model = self._get_model()

From 51a7b35a725d09d50cd3b0392901460e05726def Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Fri, 25 Sep 2020 13:46:25 -0400
Subject: [PATCH 04/14] modified whitespaces

---
 copulas/univariate/hyper_laplace.py | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index 804c93ea..565c541e 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -5,6 +5,7 @@
 from copulas.univariate.base import BoundedType, ParametricType, Univariate
 from scipy.special import gamma, psi
 
+
 def gamma_d(x):
     """
     Compute the derivative of gamma function
@@ -18,6 +19,7 @@ def gamma_d(x):
     """
     return gamma(x) * psi(x)
 
+
 def g3g1_over_g22(x):
     """
     Commute the value of Gamma(3x)Gamma(x)/Gamma^2(2x)
@@ -31,6 +33,7 @@ def g3g1_over_g22(x):
     """
     return gamma(3 * x) * gamma(x) / gamma(2 * x) ** 2
 
+
 def g3g1_over_g22_d(x):
     """
     Compute the derivative of Gamma(3x)Gamma(x)/Gamma^2(2x)
@@ -45,13 +48,14 @@ def g3g1_over_g22_d(x):
     x1 = gamma(x)
     x2 = gamma(2 * x)
     x3 = gamma(3 * x)
-    
+
     x3_d = 3 * gamma_d(3 * x) * x1 * (x2 ** 2)
     x1_d = x3 * gamma_d(x) * (x2 ** 2)
     x2_d = 4 * x3 * x1 * gamma_d(2 * x) * x2
     return (x3_d + x1_d - x2_d) / (x2 ** 4)
 
-def solve_gamma_equation(y, num_iter = 10, upper_threshold = 2, lower_threshold = 0.5):
+
+def solve_gamma_equation(y, num_iter=10, upper_threshold=2, lower_threshold=0.5):
     """
     Solve for the equation Gamma(3x)Gamma(x)/Gamma^2(2x) = y using Newton's method
 
@@ -60,7 +64,7 @@ def solve_gamma_equation(y, num_iter = 10, upper_threshold = 2, lower_threshold
 
             num_iter (int):
                 number of iterations when using Newton's method
-            
+
             upper_threshold (float):
                 Maximum value of x_new/x_old allowed in Newton's method. Should be > 1
 
@@ -89,7 +93,8 @@ class HyperLaplace(Univariate):
 
     Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html
 
-    Math derivation: HyperLaplace(k, alpha) = (gamma(loc = 0, scale = 1/k, a = 1/alpha))**(1/alpha) * Unif({-1,1})
+    Math derivation: HyperLaplace(k, alpha) = 
+    (gamma(loc = 0, scale = 1/k, a = 1/alpha))**(1/alpha) * Unif({-1,1})
     """
 
     PARAMETRIC = ParametricType.PARAMETRIC
@@ -99,7 +104,7 @@ class HyperLaplace(Univariate):
 
     _params = None
     _model = None
-    
+
     alpha = None
     k = None
 
@@ -120,7 +125,7 @@ def probability_density(self, X):
                 if the model is not fitted.
         """
         self.check_fit()
-        return self._model.pdf(abs(X) ** self.alpha) * 0.5 * alpha * abs(X) ** (self.alpha - 1)
+        return self._model.pdf(abs(X) ** self.alpha) * 0.5 * self.alpha * abs(X) ** (self.alpha - 1)
 
     def log_probability_density(self, X):
         """Compute the log of the probability density for each point in X.
@@ -194,7 +199,7 @@ def percent_point(self, U):
 
         return gamma_ppf
 
-    def sample(self, n_samples = 1):
+    def sample(self, n_samples=1):
         """Sample values from this model.
 
         Argument:
@@ -220,7 +225,8 @@ def sample(self, n_samples = 1):
     def _fit(self, X):
         """Fit the model to a non-constant random variable.
 
-        This fitting method is implemented by matching the theoritical mean and variance with the emperical mean and variance.
+        This fitting method is implemented by matching the theoritical mean 
+        and variance with the emperical mean and variance.
 
         Arguments:
             X (numpy.ndarray):
@@ -228,7 +234,7 @@ def _fit(self, X):
         """
         mean = np.mean(abs(X))
         square = np.mean(X ** 2)
-        a = solve_gamma_equation(square / (mean **2 ))
+        a = solve_gamma_equation(square / (mean ** 2))
         self.alpha = 1 / a
         self.k = (gamma(2 * a) / (gamma(a) * mean)) ** self.alpha
 

From fafaa5338300de184774708603ea9d4e8a9e4c94 Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Fri, 25 Sep 2020 13:56:02 -0400
Subject: [PATCH 05/14] modified whitespaces

---
 copulas/univariate/hyper_laplace.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index 565c541e..86611cb4 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -93,7 +93,7 @@ class HyperLaplace(Univariate):
 
     Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html
 
-    Math derivation: HyperLaplace(k, alpha) = 
+    Math derivation: HyperLaplace(k, alpha) =
     (gamma(loc = 0, scale = 1/k, a = 1/alpha))**(1/alpha) * Unif({-1,1})
     """
 
@@ -125,7 +125,7 @@ def probability_density(self, X):
                 if the model is not fitted.
         """
         self.check_fit()
-        return self._model.pdf(abs(X) ** self.alpha) * 0.5 * self.alpha * abs(X) ** (self.alpha - 1)
+        return self._model.pdf(abs(X) ** self.alpha) / 2 * self.alpha * abs(X) ** (self.alpha - 1)
 
     def log_probability_density(self, X):
         """Compute the log of the probability density for each point in X.
@@ -225,7 +225,7 @@ def sample(self, n_samples=1):
     def _fit(self, X):
         """Fit the model to a non-constant random variable.
 
-        This fitting method is implemented by matching the theoritical mean 
+        This fitting method is implemented by matching the theoritical mean
         and variance with the emperical mean and variance.
 
         Arguments:

From 0d780f71784c31a62b7a904516064bf9f85dd8b6 Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Fri, 25 Sep 2020 14:11:32 -0400
Subject: [PATCH 06/14] Regrouped imports

---
 copulas/univariate/__init__.py      | 2 +-
 copulas/univariate/hyper_laplace.py | 2 +-
 tests/unit/univariate/test_base.py  | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/copulas/univariate/__init__.py b/copulas/univariate/__init__.py
index 399bf3d3..477c86bc 100644
--- a/copulas/univariate/__init__.py
+++ b/copulas/univariate/__init__.py
@@ -3,11 +3,11 @@
 from copulas.univariate.gamma import GammaUnivariate
 from copulas.univariate.gaussian import GaussianUnivariate
 from copulas.univariate.gaussian_kde import GaussianKDE
+from copulas.univariate.hyper_laplace import HyperLaplace
 from copulas.univariate.log_laplace import LogLaplace
 from copulas.univariate.student_t import StudentTUnivariate
 from copulas.univariate.truncated_gaussian import TruncatedGaussian
 from copulas.univariate.uniform import UniformUnivariate
-from copulas.univariate.hyper_laplace import HyperLaplace
 
 __all__ = (
     'BetaUnivariate',
diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index 86611cb4..bff47eb6 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -1,9 +1,9 @@
 import numpy as np
 import scipy.stats
 from random import random
+from scipy.special import gamma, psi
 
 from copulas.univariate.base import BoundedType, ParametricType, Univariate
-from scipy.special import gamma, psi
 
 
 def gamma_d(x):
diff --git a/tests/unit/univariate/test_base.py b/tests/unit/univariate/test_base.py
index 16c2f2b7..adaaf720 100644
--- a/tests/unit/univariate/test_base.py
+++ b/tests/unit/univariate/test_base.py
@@ -7,11 +7,11 @@
 from copulas.univariate.gamma import GammaUnivariate
 from copulas.univariate.gaussian import GaussianUnivariate
 from copulas.univariate.gaussian_kde import GaussianKDE
+from copulas.univariate.hyper_laplace import HyperLaplace
 from copulas.univariate.log_laplace import LogLaplace
 from copulas.univariate.student_t import StudentTUnivariate
 from copulas.univariate.truncated_gaussian import TruncatedGaussian
 from copulas.univariate.uniform import UniformUnivariate
-from copulas.univariate.hyper_laplace import HyperLaplace
 from tests import compare_nested_iterables
 
 

From 0996c88d457de039356c2b9ff493bed23f2d368e Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Fri, 25 Sep 2020 14:57:27 -0400
Subject: [PATCH 07/14] Reordered import in hyper_laplace

---
 copulas/univariate/hyper_laplace.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index bff47eb6..fee0a958 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -1,7 +1,8 @@
-import numpy as np
-import scipy.stats
 from random import random
+
+import numpy as np
 from scipy.special import gamma, psi
+import scipy.stats
 
 from copulas.univariate.base import BoundedType, ParametricType, Univariate
 

From 5d195a2873fafa2af3e3ec3c9046a5544360d0a0 Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Fri, 25 Sep 2020 15:06:52 -0400
Subject: [PATCH 08/14] Reordered imports in hyper_laplace

---
 copulas/univariate/hyper_laplace.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index fee0a958..3e3dd2d7 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -1,8 +1,8 @@
 from random import random
 
 import numpy as np
-from scipy.special import gamma, psi
 import scipy.stats
+from scipy.special import gamma, psi
 
 from copulas.univariate.base import BoundedType, ParametricType, Univariate
 

From e40024384173d2145aa001ce22db6af70f5d79ab Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Mon, 28 Sep 2020 21:48:02 -0400
Subject: [PATCH 09/14] Implemented unit tests for hyper laplace

---
 copulas/univariate/hyper_laplace.py         |  2 +-
 tests/unit/univariate/test_hyper_laplace.py | 44 +++++++++++++++++++++
 2 files changed, 45 insertions(+), 1 deletion(-)
 create mode 100644 tests/unit/univariate/test_hyper_laplace.py

diff --git a/copulas/univariate/hyper_laplace.py b/copulas/univariate/hyper_laplace.py
index 3e3dd2d7..86431a15 100644
--- a/copulas/univariate/hyper_laplace.py
+++ b/copulas/univariate/hyper_laplace.py
@@ -23,7 +23,7 @@ def gamma_d(x):
 
 def g3g1_over_g22(x):
     """
-    Commute the value of Gamma(3x)Gamma(x)/Gamma^2(2x)
+    Compute the value of Gamma(3x)Gamma(x)/Gamma^2(2x)
 
     Arguments:
             x (float)
diff --git a/tests/unit/univariate/test_hyper_laplace.py b/tests/unit/univariate/test_hyper_laplace.py
new file mode 100644
index 00000000..6f0162cf
--- /dev/null
+++ b/tests/unit/univariate/test_hyper_laplace.py
@@ -0,0 +1,44 @@
+from unittest import TestCase
+
+import numpy as np
+from scipy.stats import norm
+
+from copulas.univariate.hyper_laplace import HyperLaplace
+
+
+class TestHyperLaplace(TestCase):
+
+    def test__fit_constant(self):
+        distribution = HyperLaplace()
+
+        distribution._fit_constant(np.array([1, 1, 1, 1]))
+
+        assert distribution._params == {
+            'loc': 1,
+            'scale': 0,
+            'a': 1
+        }
+
+    def test__fit(self): 
+        distribution = HyperLaplace()
+
+        data = norm.rvs(size=100000, loc=0, scale=1)
+        distribution._fit(data)
+
+        assert distribution._params['loc'] == 0
+        assert distribution._params['scale'] >= 1.90 and distribution._params['scale'] <= 2.10
+        assert distribution._params['a'] >= 0.48 and distribution._params['a'] <= 0.52
+
+    def test__is_constant_true(self):
+        distribution = HyperLaplace()
+
+        distribution.fit(np.array([1, 1, 1, 1]))
+
+        assert distribution._is_constant()
+
+    def test__is_constant_false(self):
+        distribution = HyperLaplace()
+
+        distribution.fit(np.array([1, 2, 3, 4]))
+
+        assert not distribution._is_constant()

From 6a6108dbb117a96e44bbf32d8920531e94e815ca Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Mon, 28 Sep 2020 21:59:44 -0400
Subject: [PATCH 10/14] deleted trailing space

---
 tests/unit/univariate/test_hyper_laplace.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/unit/univariate/test_hyper_laplace.py b/tests/unit/univariate/test_hyper_laplace.py
index 6f0162cf..7f2bf7d5 100644
--- a/tests/unit/univariate/test_hyper_laplace.py
+++ b/tests/unit/univariate/test_hyper_laplace.py
@@ -19,7 +19,7 @@ def test__fit_constant(self):
             'a': 1
         }
 
-    def test__fit(self): 
+    def test__fit(self):
         distribution = HyperLaplace()
 
         data = norm.rvs(size=100000, loc=0, scale=1)

From 670e2ccc1e8d420f17fbe62c0b853ae9e52269ee Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Mon, 28 Sep 2020 23:00:51 -0400
Subject: [PATCH 11/14] retrigger checks


From 98e13b0e4d2918e08ab8100616e5e6446f818183 Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Mon, 28 Sep 2020 23:08:57 -0400
Subject: [PATCH 12/14] retrigger checks


From 2fe1d941e3032953c726739b65363ae2432a360d Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Mon, 28 Sep 2020 23:23:17 -0400
Subject: [PATCH 13/14] retrigger checks


From 3e5c9c584bc9cc673a2666207e6683ba30baca7c Mon Sep 17 00:00:00 2001
From: Zhuofan Xie <zhuofanx@mit.edu>
Date: Mon, 28 Sep 2020 23:28:52 -0400
Subject: [PATCH 14/14] retrigger checks