Automatic Differentiation

delicatessen/estimating_equations/dose_response.py

@@ -85,6 +85,7 @@ def ee_4p_logistic(theta, X, y):
     To summarize the recommendations, be sure to examine your data (e.g., scatterplot). This will help to determine the
     initial starting values for the root-finding procedure. Otherwise, you may come across a convergence error.
 
+
     >>> estr = MEstimator(psi, init=[np.min(resp_data),
     >>>                              (np.max(resp_data)+np.min(resp_data)) / 2,
     >>>                              (np.max(resp_data)+np.min(resp_data)) / 2,
@@ -120,18 +121,19 @@ def ee_4p_logistic(theta, X, y):
     # Generalized 4PL model function for y-hat
     fx = theta[0] + (theta[3] - theta[0]) / (1 + rho)
 
-    # Using a special implementatin of natural log here
-    nested_log = np.log(X / theta[1],             # ... to avoid dose=0 issues only take log
-                        where=0 < X)              # ... where dose>0 (otherwise puts zero in place)
+    # Using a special implementation of natural log here
+    # nested_log = np.log(X / theta[1],             # ... to avoid dose=0 issues only take log
+    #                     where=0 < X)              # ... where dose>0 (otherwise puts zero in place)
+    nested_log = np.where(X > 0, np.log(X / theta[1]), 0)  # Handling when dose = 0
 
     # Calculate the derivatives for the gradient
-    deriv = np.array((1 - 1/(1 + rho),                                         # Gradient for lower limit
-                     (theta[3]-theta[0])*theta[2]/theta[1]*rho/(1+rho)**2,     # Gradient for steepness
-                     (theta[3] - theta[0]) * nested_log * rho / (1 + rho)**2,  # Gradient for ED50
-                     1 / (1 + rho)), )                                         # Gradient for upper limit
+    deriv = np.array((1 - 1/(1 + rho),                                          # Gradient for lower limit
+                     (theta[3] - theta[0])*theta[2]/theta[1]*rho/(1 + rho)**2,  # Gradient for steepness
+                     (theta[3] - theta[0])*nested_log*rho/(1 + rho)**2,         # Gradient for ED50
+                     1 / (1 + rho)), )                                          # Gradient for upper limit
 
     # Compute gradient and return for each i
-    return -2*(y-fx)*deriv
+    return -2*(y - fx)*deriv
 
 
 def ee_3p_logistic(theta, X, y, lower):
@@ -229,23 +231,8 @@ def ee_3p_logistic(theta, X, y, lower):
     Inderjit, Streibig JC, & Olofsdotter M. (2002). Joint action of phenolic acid mixtures and its significance in
     allelopathy research. *Physiologia Plantarum*, 114(3), 422-428.
     """
-    # Creating rho to cut down on typing
-    rho = (X / theta[0])**theta[1]
-
-    # Generalized 3PL model function for y-hat
-    fx = lower + (theta[2] - lower) / (1 + rho)
-
-    # Using a special implementation of natural log here
-    nested_log = np.log(X / theta[0],             # ... to avoid dose=0 issues only take log
-                        where=0 < X)              # ... where dose>0 (otherwise puts zero in place)
-
-    # Calculate the derivatives for the gradient
-    deriv = np.array(((theta[2]-lower)*theta[1]/theta[0]*rho/(1+rho)**2,     # Gradient for steepness
-                      (theta[2]-lower) * nested_log * rho / (1+rho)**2,      # Gradient for ED50
-                      1 / (1 + rho)), )                                      # Gradient for upper limit
-
-    # Compute gradient and return for each i
-    return -2*(y - fx)*deriv
+    ee_dr = ee_4p_logistic(theta=[lower, theta[0], theta[1], theta[2]], X=X, y=y)
+    return ee_dr[1:, :]
 
 
 def ee_2p_logistic(theta, X, y, lower, upper):
@@ -343,22 +330,8 @@ def ee_2p_logistic(theta, X, y, lower, upper):
     Inderjit, Streibig JC, & Olofsdotter M. (2002). Joint action of phenolic acid mixtures and its significance in
     allelopathy research. *Physiologia Plantarum*, 114(3), 422-428.
     """
-    # Creating rho to cut down on typing
-    rho = (X / theta[0])**theta[1]
-
-    # Generalized 3PL model function for y-hat
-    fx = lower + (upper - lower) / (1 + rho)
-
-    # Using a special implementatin of natural log here
-    nested_log = np.log(X / theta[0],             # ... to avoid dose=0 issues only take log
-                        where=0 < X)              # ... where dose>0 (otherwise puts zero in place)
-
-    # Calculate the derivatives for the gradient
-    deriv = np.array(((upper-lower)*theta[1]/theta[0]*rho/(1+rho)**2,     # Gradient for steepness
-                      (upper-lower) * nested_log * rho / (1+rho)**2), )   # Gradient for ED50
-
-    # Compute gradient and return for each i
-    return -2*(y-fx)*deriv
+    ee_dr = ee_4p_logistic(theta=[lower, theta[0], theta[1], upper], X=X, y=y)
+    return ee_dr[1:3, :]
 
 
 def ee_effective_dose_delta(theta, y, delta, steepness, ed50, lower, upper):

delicatessen/estimating_equations/regression.py

@@ -1,11 +1,11 @@
 import warnings
 import numpy as np
 from scipy.stats import norm, cauchy
-from scipy.special import digamma, polygamma
 
 from delicatessen.utilities import (logit, inverse_logit, identity,
                                     robust_loss_functions,
-                                    additive_design_matrix)
+                                    additive_design_matrix,
+                                    digamma, polygamma, standard_normal_cdf, standard_normal_pdf)
 from delicatessen.estimating_equations.processing import generate_weights
 
 #################################################################
@@ -76,7 +76,7 @@ def ee_regression(theta, X, y, model, weights=None, offset=None):
     >>> data['Z'] = np.random.normal(size=n)
     >>> data['Y1'] = 0.5 + 2*data['X'] - 1*data['Z'] + np.random.normal(loc=0, size=n)
     >>> data['Y2'] = np.random.binomial(n=1, p=logistic.cdf(0.5 + 2*data['X'] - 1*data['Z']), size=n)
-    >>> data['Y3'] = np.random.poisson(np.exp(lam=0.5 + 2*data['X'] - 1*data['Z']), size=n)
+    >>> data['Y3'] = np.random.poisson(lam=np.exp(0.5 + 2*data['X'] - 1*data['Z']), size=n)
     >>> data['C'] = 1
 
     Note that ``C`` here is set to all 1's. This will be the intercept in the regression.
@@ -1351,11 +1351,15 @@ def _inverse_link_(betax, link):
         py = 1 - np.exp(-1*np.exp(betax))       # Inverse link
         dpy = np.exp(betax - np.exp(betax))     # Derivative of inverse link
     elif link == 'probit':
-        py = norm.cdf(betax)                    # Inverse link
-        dpy = norm.pdf(betax)                   # Derivative of inverse link
+        # py = norm.cdf(betax)                  # Inverse link
+        # dpy = norm.pdf(betax)                 # Derivative of inverse link
+        py = standard_normal_cdf(x=betax)       # Inverse link
+        dpy = standard_normal_pdf(x=betax)      # Derivative of inverse link
     elif link in ['cauchit', 'cauchy']:
-        py = cauchy.cdf(betax)                  # Inverse link
-        dpy = cauchy.pdf(betax)                 # Derivative of inverse link
+        # py = cauchy.cdf(betax)                # Inverse link
+        # dpy = cauchy.pdf(betax)               # Derivative of inverse link
+        py = (1/np.pi)*np.arctan(betax) + 0.5   # Inverse link
+        dpy = 1 / (np.pi*(1 + betax**2))        # Derivative of inverse link (by-hand)
     elif link in ['square_root', 'sqrt']:
         py = betax**2                           # Inverse link
         dpy = 2 * betax                         # Derivative of inverse link