Skip to content
This repository has been archived by the owner on Oct 31, 2024. It is now read-only.

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history
for more information, see https://pre-commit.ci
  • Loading branch information
pre-commit-ci[bot] committed Oct 30, 2023
1 parent 3ac554b commit e467148
Show file tree
Hide file tree
Showing 15 changed files with 61 additions and 85 deletions.
92 changes: 46 additions & 46 deletions development/hessian_approximation/compute_hessian.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 = (
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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
)
)
Expand All @@ -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
)
)
Expand Down
4 changes: 0 additions & 4 deletions development/tests/property/property_auxiliary.py
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,6 @@ def print_rslt_ext(start, timeout, rslt, err_msg):
current_time = datetime.now().strftime("%Y-%m-%d %H:%M:%S")

with open("property.grmpy.info", "w") as outfile:

# Write out some header information.
outfile.write("\n\n")
str_ = "\t{0[0]:<15}{0[1]:<20}\n\n"
Expand All @@ -64,14 +63,12 @@ def print_rslt_ext(start, timeout, rslt, err_msg):
modules = sorted(rslt.keys())

for module in modules:

outfile.write("\n " + module.replace(".py", "") + "\n\n")

string = "{:>18}{:>15}{:>15}\n\n"
outfile.write(string.format("Test", "Success", "Failure"))

for test in sorted(rslt[module].keys()):

stat = rslt[module][test]

string = "{:>18}{:>15}{:>15}\n"
Expand All @@ -81,7 +78,6 @@ def print_rslt_ext(start, timeout, rslt, err_msg):
outfile.write("\n" + "-" * 79 + "\n\n")

for err in err_msg:

module, test, seed, msg = err

string = "MODULE {:<25} TEST {:<15} SEED {:<15}\n\n"
Expand Down
2 changes: 0 additions & 2 deletions development/tests/property/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,6 @@ def run(args):
print_rslt(rslt, err_msg)

while True:

seed = random.randrange(1, 100000)
dirname = get_random_string()
np.random.seed(seed)
Expand Down Expand Up @@ -86,7 +85,6 @@ def run(args):


if __name__ == "__main__":

args = process_command_line_arguments("property")

run(args)
1 change: 0 additions & 1 deletion development/tests/regression/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,6 @@ def check_vault(num_tests=100):


if __name__ == "__main__":

parser = argparse.ArgumentParser(
description="Work with regression tests for package."
)
Expand Down
2 changes: 0 additions & 2 deletions development/tests/reliability/reliability.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,6 @@ def monte_carlo(file, grid_points):

# Loop over different correlations between V and U_1
for rho in np.linspace(0.00, 0.99, grid_points):

# Readjust the initialization file values to add correlation
model_spec = read(file)
sim_spec = read("reliability.grmpy.yml")
Expand Down Expand Up @@ -153,7 +152,6 @@ def create_plots(effects, true):


if __name__ == "__main__":

ATE = get_effect_grmpy("reliability.grmpy.yml")

x = monte_carlo("reliability.grmpy.yml", 10)
Expand Down
1 change: 0 additions & 1 deletion development/tests/replication/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,6 @@ def calculate_cof_int(rslt, init_dict, data_frame, mte, quantiles):


if __name__ == "__main__":

init_dict = read("replication.grmpy.yml")
# Estimate the coefficients
rslt = fit("replication.grmpy.yml")
Expand Down
1 change: 0 additions & 1 deletion docs/source/figures/create.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,5 +21,4 @@ def create_figures():


if __name__ == "__main__":

create_figures()
1 change: 0 additions & 1 deletion docs/source/figures/scripts/fig-reliability-par.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,5 @@ def plot_rslts(rslt, file):


if __name__ == "__main__":

rslt_dict = fit(RESOURCE_DIR + "/replication.grmpy.yml")
plot_rslts(rslt_dict, RESOURCE_DIR + "/replication.grmpy.yml")
1 change: 0 additions & 1 deletion docs/source/figures/scripts/fig-treatment-effects.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@


def plot_treatment_effect():

x_axis = np.arange(-2, 4, 0.001)
ax = plt.figure(figsize=(14, 6))

Expand Down
Loading

0 comments on commit e467148

Please sign in to comment.