[pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

Author: pre-commit-ci[bot]

development/hessian_approximation/compute_hessian.py

@@ -29,9 +29,9 @@ def compute_hessian(x0, X1, X0, Z1, Z0, Y1, Y0):
 
     # aux_params
     nu1 = (Y1 - np.dot(beta1, X1.T)) / sd1
-    lambda1 = (np.dot(gamma, Z1.T) - rho1v * nu1) / (np.sqrt(1 - rho1v ** 2))
+    lambda1 = (np.dot(gamma, Z1.T) - rho1v * nu1) / (np.sqrt(1 - rho1v**2))
     nu0 = (Y0 - np.dot(beta0, X0.T)) / sd0
-    lambda0 = (np.dot(gamma, Z0.T) - rho0v * nu0) / (np.sqrt(1 - rho0v ** 2))
+    lambda0 = (np.dot(gamma, Z0.T) - rho0v * nu0) / (np.sqrt(1 - rho0v**2))
 
     eta1 = (
         -lambda1 * norm.pdf(lambda1) * norm.cdf(lambda1) - norm.pdf(lambda1) ** 2
@@ -113,21 +113,21 @@ def calc_hess_beta1(
     """
 
     # define some auxiliary variables
-    rho_aux1 = lambda1 * rho1v / (1 - rho1v ** 2) - nu1 / (1 - rho1v ** 2) ** 0.5
-    rho_aux2 = rho1v ** 2 / ((1 - rho1v ** 2) ** (3 / 2)) + 1 / (1 - rho1v ** 2) ** 0.5
-    sd_aux1 = rho1v ** 2 / (1 - rho1v ** 2)
-    sd_aux2 = rho1v / np.sqrt(1 - rho1v ** 2)
+    rho_aux1 = lambda1 * rho1v / (1 - rho1v**2) - nu1 / (1 - rho1v**2) ** 0.5
+    rho_aux2 = rho1v**2 / ((1 - rho1v**2) ** (3 / 2)) + 1 / (1 - rho1v**2) ** 0.5
+    sd_aux1 = rho1v**2 / (1 - rho1v**2)
+    sd_aux2 = rho1v / np.sqrt(1 - rho1v**2)
     # derivation wrt beta1
     der_b1_beta1 = -(
-        X1X1 * (rho1v ** 2 / (1 - rho1v ** 2)) * 1 / sd1 ** 2 - X1.T @ X1 / sd1 ** 2
+        X1X1 * (rho1v**2 / (1 - rho1v**2)) * 1 / sd1**2 - X1.T @ X1 / sd1**2
     )
     # add zeros for derv beta 0
     der_b1_beta1 = np.concatenate(
         (der_b1_beta1, np.zeros((n_col_X1, n_col_X0))), axis=1
     )
 
     # derivation wrt gamma
-    der_b1_gamma = -(X1Z1 * rho1v / (sd1 * (1 - rho1v ** 2)))
+    der_b1_gamma = -(X1Z1 * rho1v / (sd1 * (1 - rho1v**2)))
     der_b1_gamma = np.concatenate((der_b1_beta1, der_b1_gamma), axis=1)
     # derv wrt sigma 1
     der_b1_sd = (
@@ -150,7 +150,7 @@ def calc_hess_beta1(
     # derv wrt rho1
     der_b1_rho = (
         -(
-            eta1 * rho_aux1 * rho1v / ((1 - rho1v ** 2) ** 0.5)
+            eta1 * rho_aux1 * rho1v / ((1 - rho1v**2) ** 0.5)
             + norm.pdf(lambda1) / norm.cdf(lambda1) * rho_aux2
         )
         * 1
@@ -176,21 +176,21 @@ def calc_hess_beta0(
     """
 
     # define some aux_vars
-    rho_aux1 = lambda0 * rho0v / (1 - rho0v ** 2) - nu0 / (1 - rho0v ** 2) ** 0.5
-    rho_aux2 = rho0v ** 2 / ((1 - rho0v ** 2) ** (3 / 2)) + 1 / (1 - rho0v ** 2) ** 0.5
-    sd_aux1 = rho0v ** 2 / (1 - rho0v ** 2)
-    sd_aux2 = rho0v / (np.sqrt(1 - rho0v ** 2))
+    rho_aux1 = lambda0 * rho0v / (1 - rho0v**2) - nu0 / (1 - rho0v**2) ** 0.5
+    rho_aux2 = rho0v**2 / ((1 - rho0v**2) ** (3 / 2)) + 1 / (1 - rho0v**2) ** 0.5
+    sd_aux1 = rho0v**2 / (1 - rho0v**2)
+    sd_aux2 = rho0v / (np.sqrt(1 - rho0v**2))
 
     # add zeros for beta0
     der_b0_beta1 = np.zeros((n_col_X1, n_col_X0))
 
     # beta0
     der_b0_beta0 = (
-        -(X0X0 * (rho0v ** 2 / (1 - rho0v ** 2)) * 1 / sd0 ** 2) + X0.T @ X0 / sd0 ** 2
+        -(X0X0 * (rho0v**2 / (1 - rho0v**2)) * 1 / sd0**2) + X0.T @ X0 / sd0**2
     )
     der_b0_beta0 = np.concatenate((der_b0_beta1, der_b0_beta0), axis=1)
     # gamma
-    der_b0_gamma = -X0Z0 * rho0v / (1 - rho0v ** 2) * 1 / sd0
+    der_b0_gamma = -X0Z0 * rho0v / (1 - rho0v**2) * 1 / sd0
     der_b0_gamma = np.concatenate((der_b0_beta0, der_b0_gamma), axis=1)
 
     # add zeros for sigma1 and rho1
@@ -204,15 +204,15 @@ def calc_hess_beta0(
             - 2 * nu0
         )
         * 1
-        / sd0 ** 2
+        / sd0**2
     )
     der_b0_sd = np.expand_dims((der_b0_sd.T @ X0), 1)
     der_b0_sd = np.concatenate((der_b0_gamma, der_b0_sd), axis=1)
 
     # rho
     der_b0_rho = (
         (
-            eta0 * -rho_aux1 * (rho0v / ((1 - rho0v ** 2) ** 0.5))
+            eta0 * -rho_aux1 * (rho0v / ((1 - rho0v**2) ** 0.5))
             + norm.pdf(lambda0) / (1 - norm.cdf(lambda0)) * rho_aux2
         )
         * 1
@@ -252,7 +252,7 @@ def calc_hess_gamma(
 
     der_gamma_beta = np.zeros((Z1.shape[1], num_col_X1X0))
 
-    der_g_gamma = -(1 / (1 - rho1v ** 2) * Z1Z1 + 1 / (1 - rho0v ** 2) * Z0Z0)
+    der_g_gamma = -(1 / (1 - rho1v**2) * Z1Z1 + 1 / (1 - rho0v**2) * Z0Z0)
     der_g_gamma = np.concatenate((der_gamma_beta, der_g_gamma), axis=1)
 
     # sigma1
@@ -261,19 +261,19 @@ def calc_hess_gamma(
         @ np.einsum("ij, i ->ij", X1, nu1)
         / sd1
         * rho1v
-        / (1 - rho1v ** 2)
+        / (1 - rho1v**2)
     )[:, 0]
     der_g_sd1 = np.expand_dims(der_g_sd1, 0).T
     der_g_sd1 = np.concatenate((der_g_gamma, der_g_sd1), axis=1)
 
     # rho1
     aux_rho11 = np.einsum("ij, i ->ij", Z1, eta1).T @ (
-        lambda1 * rho1v / (1 - rho1v ** 2) - nu1 / np.sqrt(1 - rho1v ** 2)
+        lambda1 * rho1v / (1 - rho1v**2) - nu1 / np.sqrt(1 - rho1v**2)
     )
     aux_rho21 = Z1.T @ (norm.pdf(lambda1) / norm.cdf(lambda1))
 
-    der_g_rho1 = -aux_rho11 * 1 / (np.sqrt(1 - rho1v ** 2)) - aux_rho21 * rho1v / (
-        (1 - rho1v ** 2) ** (3 / 2)
+    der_g_rho1 = -aux_rho11 * 1 / (np.sqrt(1 - rho1v**2)) - aux_rho21 * rho1v / (
+        (1 - rho1v**2) ** (3 / 2)
     )
     der_g_rho1 = np.expand_dims(der_g_rho1, 0).T
     der_g_rho1 = np.concatenate((der_g_sd1, der_g_rho1), axis=1)
@@ -284,19 +284,19 @@ def calc_hess_gamma(
         @ np.einsum("ij, i ->ij", X0, nu0)
         / sd0
         * rho0v
-        / (1 - rho0v ** 2)
+        / (1 - rho0v**2)
     )[:, 0]
     der_g_sd0 = np.expand_dims(der_g_sd0, 0).T
     der_g_sd0 = np.concatenate((der_g_rho1, -der_g_sd0), axis=1)
 
     # rho1
     aux_rho10 = np.einsum("ij, i ->ij", Z0, eta0).T @ (
-        lambda0 * rho0v / (1 - rho0v ** 2) - nu0 / np.sqrt(1 - rho0v ** 2)
+        lambda0 * rho0v / (1 - rho0v**2) - nu0 / np.sqrt(1 - rho0v**2)
     )
     aux_rho20 = -Z0.T @ (norm.pdf(lambda0) / (1 - norm.cdf(lambda0)))
 
-    der_g_rho0 = aux_rho10 * 1 / (np.sqrt(1 - rho0v ** 2)) + aux_rho20 * rho0v / (
-        (1 - rho0v ** 2) ** (3 / 2)
+    der_g_rho0 = aux_rho10 * 1 / (np.sqrt(1 - rho0v**2)) + aux_rho20 * rho0v / (
+        (1 - rho0v**2) ** (3 / 2)
     )
     der_g_rho0 = np.expand_dims(-der_g_rho0, 0).T
     der_g_rho0 = np.concatenate((der_g_sd0, der_g_rho0), axis=1)
@@ -328,52 +328,52 @@ def calc_hess_dist(
     Delta_sd1 = (
         +1 / sd1
         - (norm.pdf(lambda1) / norm.cdf(lambda1))
-        * (rho1v * nu1 / (np.sqrt(1 - rho1v ** 2) * sd1))
-        - nu1 ** 2 / sd1
+        * (rho1v * nu1 / (np.sqrt(1 - rho1v**2) * sd1))
+        - nu1**2 / sd1
     )
     Delta_sd1_der = (
         nu1
         / sd1
         * (
-            -eta1 * (rho1v ** 2 * nu1) / (1 - rho1v ** 2)
-            + (norm.pdf(lambda1) / norm.cdf(lambda1)) * rho1v / np.sqrt(1 - rho1v ** 2)
+            -eta1 * (rho1v**2 * nu1) / (1 - rho1v**2)
+            + (norm.pdf(lambda1) / norm.cdf(lambda1)) * rho1v / np.sqrt(1 - rho1v**2)
             + 2 * nu1
         )
     )
 
     Delta_sd0 = (
         +1 / sd0
         + (norm.pdf(lambda0) / (1 - norm.cdf(lambda0)))
-        * (rho0v * nu0 / (np.sqrt(1 - rho0v ** 2) * sd0))
-        - nu0 ** 2 / sd0
+        * (rho0v * nu0 / (np.sqrt(1 - rho0v**2) * sd0))
+        - nu0**2 / sd0
     )
     Delta_sd0_der = (
         nu0
         / sd0
         * (
-            -eta0 * (rho0v ** 2 * nu0) / (1 - rho0v ** 2)
+            -eta0 * (rho0v**2 * nu0) / (1 - rho0v**2)
             - (norm.pdf(lambda0) / (1 - norm.cdf(lambda0)))
             * rho0v
-            / np.sqrt(1 - rho0v ** 2)
+            / np.sqrt(1 - rho0v**2)
             + 2 * nu0
         )
     )
 
-    aux_rho11 = lambda1 * rho1v / (1 - rho1v ** 2) - nu1 / np.sqrt(1 - rho1v ** 2)
-    aux_rho12 = 1 / (1 - rho1v ** 2) ** (3 / 2)
+    aux_rho11 = lambda1 * rho1v / (1 - rho1v**2) - nu1 / np.sqrt(1 - rho1v**2)
+    aux_rho12 = 1 / (1 - rho1v**2) ** (3 / 2)
 
-    aux_rho_rho11 = (np.dot(gamma, Z1.T) * rho1v - nu1) / (1 - rho1v ** 2) ** (3 / 2)
+    aux_rho_rho11 = (np.dot(gamma, Z1.T) * rho1v - nu1) / (1 - rho1v**2) ** (3 / 2)
     aux_rho_rho12 = (
-        2 * np.dot(gamma, Z1.T) * rho1v ** 2 + np.dot(gamma, Z1.T) - 3 * nu1 * rho1v
-    ) / (1 - rho1v ** 2) ** (5 / 2)
+        2 * np.dot(gamma, Z1.T) * rho1v**2 + np.dot(gamma, Z1.T) - 3 * nu1 * rho1v
+    ) / (1 - rho1v**2) ** (5 / 2)
 
-    aux_rho01 = lambda0 * rho0v / (1 - rho0v ** 2) - nu0 / np.sqrt(1 - rho0v ** 2)
-    aux_rho02 = 1 / (1 - rho0v ** 2) ** (3 / 2)
+    aux_rho01 = lambda0 * rho0v / (1 - rho0v**2) - nu0 / np.sqrt(1 - rho0v**2)
+    aux_rho02 = 1 / (1 - rho0v**2) ** (3 / 2)
 
-    aux_rho_rho01 = (np.dot(gamma, Z0.T) * rho0v - nu0) / (1 - rho0v ** 2) ** (3 / 2)
+    aux_rho_rho01 = (np.dot(gamma, Z0.T) * rho0v - nu0) / (1 - rho0v**2) ** (3 / 2)
     aux_rho_rho02 = (
-        2 * np.dot(gamma, Z0.T) * rho0v ** 2 + np.dot(gamma, Z0.T) - 3 * nu0 * rho0v
-    ) / (1 - rho0v ** 2) ** (5 / 2)
+        2 * np.dot(gamma, Z0.T) * rho0v**2 + np.dot(gamma, Z0.T) - 3 * nu0 * rho0v
+    ) / (1 - rho0v**2) ** (5 / 2)
 
     # for sigma1
 
@@ -385,7 +385,7 @@ def calc_hess_dist(
         1
         / sd1
         * (
-            -eta1 * aux_rho11 * (rho1v * nu1) / (np.sqrt(1 - rho1v ** 2))
+            -eta1 * aux_rho11 * (rho1v * nu1) / (np.sqrt(1 - rho1v**2))
             - (norm.pdf(lambda1) / norm.cdf(lambda1)) * aux_rho12 * nu1
         )
     )
@@ -411,7 +411,7 @@ def calc_hess_dist(
         1
         / sd0
         * (
-            -eta0 * aux_rho01 * (rho0v * nu0) / (np.sqrt(1 - rho0v ** 2))
+            -eta0 * aux_rho01 * (rho0v * nu0) / (np.sqrt(1 - rho0v**2))
             + (norm.pdf(lambda0) / (1 - norm.cdf(lambda0))) * aux_rho02 * nu0
         )
     )

grmpy/estimate/estimate_par.py

@@ -474,10 +474,10 @@ def log_likelihood(
     sd1, sd0, rho1v, rho0v = x0[-4], x0[-2], x0[-3], x0[-1]
 
     nu1 = (Y1 - np.dot(beta1, X1.T)) / sd1
-    lambda1 = (np.dot(gamma, Z1.T) - rho1v * nu1) / (np.sqrt(1 - rho1v ** 2))
+    lambda1 = (np.dot(gamma, Z1.T) - rho1v * nu1) / (np.sqrt(1 - rho1v**2))
 
     nu0 = (Y0 - np.dot(beta0, X0.T)) / sd0
-    lambda0 = (np.dot(gamma, Z0.T) - rho0v * nu0) / (np.sqrt(1 - rho0v ** 2))
+    lambda0 = (np.dot(gamma, Z0.T) - rho0v * nu0) / (np.sqrt(1 - rho0v**2))
 
     treated = (1 / sd1) * norm.pdf(nu1) * norm.cdf(lambda1)
     untreated = (1 / sd0) * norm.pdf(nu0) * (1 - norm.cdf(lambda0))
@@ -916,7 +916,7 @@ def gradient(X1, X0, Z1, Z0, nu1, nu0, lambda1, lambda0, gamma, sd1, sd0, rho1v,
             "ij, i ->ij",
             X1,
             -(norm.pdf(lambda1) / norm.cdf(lambda1))
-            * (rho1v / (np.sqrt(1 - rho1v ** 2) * sd1))
+            * (rho1v / (np.sqrt(1 - rho1v**2) * sd1))
             - nu1 / sd1,
         ),
         0,
@@ -929,7 +929,7 @@ def gradient(X1, X0, Z1, Z0, nu1, nu0, lambda1, lambda0, gamma, sd1, sd0, rho1v,
             X0,
             norm.pdf(lambda0)
             / (1 - norm.cdf(lambda0))
-            * (rho0v / (np.sqrt(1 - rho0v ** 2) * sd0))
+            * (rho0v / (np.sqrt(1 - rho0v**2) * sd0))
             - nu0 / sd0,
         ),
         0,
@@ -939,8 +939,8 @@ def gradient(X1, X0, Z1, Z0, nu1, nu0, lambda1, lambda0, gamma, sd1, sd0, rho1v,
         * (
             +1 / sd1
             - (norm.pdf(lambda1) / norm.cdf(lambda1))
-            * (rho1v * nu1 / (np.sqrt(1 - rho1v ** 2) * sd1))
-            - nu1 ** 2 / sd1
+            * (rho1v * nu1 / (np.sqrt(1 - rho1v**2) * sd1))
+            - nu1**2 / sd1
         ),
         keepdims=True,
     )
@@ -949,16 +949,16 @@ def gradient(X1, X0, Z1, Z0, nu1, nu0, lambda1, lambda0, gamma, sd1, sd0, rho1v,
         * (
             +1 / sd0
             + (norm.pdf(lambda0) / (1 - norm.cdf(lambda0)))
-            * (rho0v * nu0 / (np.sqrt(1 - rho0v ** 2) * sd0))
-            - nu0 ** 2 / sd0
+            * (rho0v * nu0 / (np.sqrt(1 - rho0v**2) * sd0))
+            - nu0**2 / sd0
         ),
         keepdims=True,
     )
     grad_rho1v = np.sum(
         (
             -(norm.pdf(lambda1) / norm.cdf(lambda1))
             * ((np.dot(gamma, Z1.T) * rho1v) - nu1)
-            / (1 - rho1v ** 2) ** (1 / 2)
+            / (1 - rho1v**2) ** (1 / 2)
         ),
         keepdims=True,
     )
@@ -967,7 +967,7 @@ def gradient(X1, X0, Z1, Z0, nu1, nu0, lambda1, lambda0, gamma, sd1, sd0, rho1v,
         (
             (norm.pdf(lambda0) / (1 - norm.cdf(lambda0)))
             * ((np.dot(gamma, Z0.T) * rho0v) - nu0)
-            / (1 - rho0v ** 2) ** (1 / 2)
+            / (1 - rho0v**2) ** (1 / 2)
         ),
         keepdims=True,
     )
@@ -976,14 +976,14 @@ def gradient(X1, X0, Z1, Z0, nu1, nu0, lambda1, lambda0, gamma, sd1, sd0, rho1v,
         np.einsum(
             "ij, i ->ij",
             Z1,
-            (norm.pdf(lambda1) / norm.cdf(lambda1)) * 1 / np.sqrt(1 - rho1v ** 2),
+            (norm.pdf(lambda1) / norm.cdf(lambda1)) * 1 / np.sqrt(1 - rho1v**2),
         )
     ) - sum(
         np.einsum(
             "ij, i ->ij",
             Z0,
             (norm.pdf(lambda0) / (1 - norm.cdf(lambda0)))
-            * (1 / np.sqrt(1 - rho0v ** 2)),
+            * (1 / np.sqrt(1 - rho0v**2)),
         )
     )
 
@@ -1042,10 +1042,10 @@ def gradient_hessian(x0, X1, X0, Z1, Z0, Y1, Y0):
     # compute gradient for beta 1
 
     nu1 = (Y1 - np.dot(beta1, X1.T)) / sd1
-    lambda1 = (np.dot(gamma, Z1.T) - rho1v * nu1) / (np.sqrt(1 - rho1v ** 2))
+    lambda1 = (np.dot(gamma, Z1.T) - rho1v * nu1) / (np.sqrt(1 - rho1v**2))
 
     nu0 = (Y0 - np.dot(beta0, X0.T)) / sd0
-    lambda0 = (np.dot(gamma, Z0.T) - rho0v * nu0) / (np.sqrt(1 - rho0v ** 2))
+    lambda0 = (np.dot(gamma, Z0.T) - rho0v * nu0) / (np.sqrt(1 - rho0v**2))
 
     grad = gradient(
         X1, X0, Z1, Z0, nu1, nu0, lambda1, lambda0, gamma, sd1, sd0, rho1v, rho0v
@@ -1054,7 +1054,7 @@ def gradient_hessian(x0, X1, X0, Z1, Z0, Y1, Y0):
     multiplier = np.concatenate(
         (
             np.ones(len(grad[:-4])),
-            np.array([1 / sd1, 1 / (1 - rho1v ** 2), 1 / sd0, 1 / (1 - rho0v ** 2)]),
+            np.array([1 / sd1, 1 / (1 - rho1v**2), 1 / sd0, 1 / (1 - rho0v**2)]),
         )
     )