First and second derivatives (wrt x_t) of log-density of Y_t|X_t=xt
195085
def diff_ft(self, xt, yt):
        
        a, b = self.a, self.b 
        ex = np.exp(a + np.matmul(b, xt))  # shape = (dy,)
        grad = (-np.sum(ex[:, np.newaxis] * b, axis=0)
                + np.sum(yt.flatten()[:, np.newaxis] * b, axis=0)) #TODO flatten
        hess = np.zeros((self.dx, self.dx))
        for k in range(self.dy):
            hess -= ex[k] * np.outer(b[k,:], b[k,:]) 
        return grad, hess
167702
def _get_gradient_log_pdf(self):
        
        sub_vec = self.variable_assignments - self.model.mean.flatten()
        grad = - np.dot(self.model.precision_matrix, sub_vec)
        log_pdf = 0.5 * np.dot(sub_vec, grad)

        return grad, log_pdf
67032
def _param_grad_helper(self, dL_dK, X, X2, target):
        """"""
        if X2 is None: X2 = X
        dist = np.abs(X - X2.T)

        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU

        dk1theta1 = np.exp(-ly*dist)*2*(-lu)/(lu+ly)**3
        #c=np.sqrt(3)
        #t1=c/lu
        #t2=1/ly
        #dk1theta1=np.exp(-dist*ly)*t2*( (2*c*t2+2*t1)/(c*t2+t1)**2 -2*(2*c*t2*t1+t1**2)/(c*t2+t1)**3   )
        
        dk2theta1 = 1*( 
            np.exp(-lu*dist)*dist*(-ly+2*lu-lu*ly*dist+dist*lu**2)*(ly-lu)**(-2) + np.exp(-lu*dist)*(-2+ly*dist-2*dist*lu)*(ly-lu)**(-2) 
            +np.exp(-dist*lu)*(ly-2*lu+ly*lu*dist-dist*lu**2)*2*(ly-lu)**(-3) 
            +np.exp(-dist*ly)*2*(ly-lu)**(-2)
            +np.exp(-dist*ly)*2*(2*lu-ly)*(ly-lu)**(-3)
            )
      
        dk3theta1 = np.exp(-dist*lu)*(lu+ly)**(-2)*((2*lu+ly+dist*lu**2+lu*ly*dist)*(-dist-2/(lu+ly))+2+2*lu*dist+ly*dist)

        dktheta1 = self.varianceU*self.varianceY*(dk1theta1+dk2theta1+dk3theta1)




        dk1theta2 = np.exp(-ly*dist) * ((lu+ly)**(-2)) * (  (-dist)*(2*lu+ly)  +  1  +  (-2)*(2*lu+ly)/(lu+ly)  )

        dk2theta2 = 1*(
            np.exp(-dist*lu)*(ly-lu)**(-2) * ( 1+lu*dist+(-2)*(ly-2*lu+lu*ly*dist-dist*lu**2)*(ly-lu)**(-1) )
            +np.exp(-dist*ly)*(ly-lu)**(-2) * ( (-dist)*(2*lu-ly) -1+(2*lu-ly)*(-2)*(ly-lu)**(-1) )
            )

        dk3theta2 = np.exp(-dist*lu) * (-3*lu-ly-dist*lu**2-lu*ly*dist)/(lu+ly)**3

        dktheta2 = self.varianceU*self.varianceY*(dk1theta2 + dk2theta2 +dk3theta2)



        k1 = np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
        k2 = (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
        k3 = np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
        dkdvar = k1+k2+k3

        target[0] += np.sum(self.varianceY*dkdvar * dL_dK)
        target[1] += np.sum(self.varianceU*dkdvar * dL_dK)
        target[2] += np.sum(dktheta1*(-np.sqrt(3)*self.lengthscaleU**(-2)) * dL_dK)
        target[3] += np.sum(dktheta2*(-self.lengthscaleY**(-2)) * dL_dK)
399652
def derivative(func, x0, dx=1.0, n=1, args=(), order=3):
    '''
    '''
    if order < n + 1:
        raise ValueError
    if order % 2 == 0:
        raise ValueError
    weights = central_diff_weights(order, n)
    tot = 0.0
    ho = order >> 1
    for k in range(order):
        tot += weights[k]*func(x0 + (k - ho)*dx, *args)
    return tot/product([dx]*n)
1088763
def derivative(self, x):
        """"""
        if self.is_linear:
            return self
        else:
            left = self.right(x) * self.left.derivative(x)
            right = self.left(x) * self.right.derivative(x)
            return left + right
1109270
def derivatives(self, x, y, theta_E, r_trunc, center_x=0, center_y=0):
        
        x_shift = x - center_x
        y_shift = y - center_y

        dphi_dr = self._dphi_dr(x_shift, y_shift, theta_E, r_trunc)
        dr_dx, dr_dy = self._dr_dx(x_shift, y_shift)
        f_x = dphi_dr * dr_dx
        f_y = dphi_dr * dr_dy
        return f_x, f_y
151891
def partial_derivative_scalar(self, U, V, y=0):
        """"""

        self.check_fit()

        X = np.column_stack((U, V))
        return self.partial_derivative(X, y)
1097762
def log_conditional_likelihood(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        ll = 0
        for t in range(self.T):
            ll += self.local_log_likelihood(x[t], self.data[t], self.inputs[t])
        return ll
588263
def density_2d(self, x, y, rho0, Rs, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        X = r/Rs
        sigma0 = self.rho2sigma(rho0, Rs)
        if isinstance(X, int) or isinstance(X, float):
            if X == 1:
                X = 1.000001
        else:
            X[X == 1] = 1.000001
        sigma = sigma0 / (X**2-1)**2 * (-3 + (2+X**2)*self._F(X))
        return sigma
324454
def black_box_variational_inference(logprob, D, num_samples):
    """"""

    def unpack_params(params):
        # Variational dist is a diagonal Gaussian.
        mean, log_std = params[:D], params[D:]
        return mean, log_std

    def gaussian_entropy(log_std):
        return 0.5 * D * (1.0 + np.log(2*np.pi)) + np.sum(log_std)

    rs = npr.RandomState(0)
    def variational_objective(params, t):
        """Provides a stochastic estimate of the variational lower bound."""
        mean, log_std = unpack_params(params)
        samples = rs.randn(num_samples, D) * np.exp(log_std) + mean
        lower_bound = gaussian_entropy(log_std) + np.mean(logprob(samples, t))
        return -lower_bound

    gradient = grad(variational_objective)

    return variational_objective, gradient, unpack_params
924908
def _ndtri(y):
    
    # approximation for 0 <= abs(z - 0.5) <= 3/8
    P0 = [
        -5.99633501014107895267E1,
        9.80010754185999661536E1,
        -5.66762857469070293439E1,
        1.39312609387279679503E1,
        -1.23916583867381258016E0,
    ]

    Q0 = [
        1.95448858338141759834E0,
        4.67627912898881538453E0,
        8.63602421390890590575E1,
        -2.25462687854119370527E2,
        2.00260212380060660359E2,
        -8.20372256168333339912E1,
        1.59056225126211695515E1,
        -1.18331621121330003142E0,
    ]

    # Approximation for interval z = sqrt(-2 log y ) between 2 and 8
    # i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
    P1 = [
        4.05544892305962419923E0,
        3.15251094599893866154E1,
        5.71628192246421288162E1,
        4.40805073893200834700E1,
        1.46849561928858024014E1,
        2.18663306850790267539E0,
        -1.40256079171354495875E-1,
        -3.50424626827848203418E-2,
        -8.57456785154685413611E-4,
    ]

    Q1 = [
        1.57799883256466749731E1,
        4.53907635128879210584E1,
        4.13172038254672030440E1,
        1.50425385692907503408E1,
        2.50464946208309415979E0,
        -1.42182922854787788574E-1,
        -3.80806407691578277194E-2,
        -9.33259480895457427372E-4,
    ]

    # Approximation for interval z = sqrt(-2 log y ) between 8 and 64
    # i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
    P2 = [
        3.23774891776946035970E0,
        6.91522889068984211695E0,
        3.93881025292474443415E0,
        1.33303460815807542389E0,
        2.01485389549179081538E-1,
        1.23716634817820021358E-2,
        3.01581553508235416007E-4,
        2.65806974686737550832E-6,
        6.23974539184983293730E-9,
    ]

    Q2 = [
        6.02427039364742014255E0,
        3.67983563856160859403E0,
        1.37702099489081330271E0,
        2.16236993594496635890E-1,
        1.34204006088543189037E-2,
        3.28014464682127739104E-4,
        2.89247864745380683936E-6,
        6.79019408009981274425E-9,
    ]

    sign_flag = 1

    if y > (1 - EXP_NEG2):
        y = 1 - y
        sign_flag = 0

    # Shortcut case where we don't need high precision
    # between -0.135 and 0.135
    if y > EXP_NEG2:
        y -= 0.5
        y2 = y ** 2
        x = y + y * (y2 * _polevl(y2, P0, 4) / _p1evl(y2, Q0, 8))
        x = x * ROOT_2PI
        return x

    x = math.sqrt(-2.0 * math.log(y))
    x0 = x - math.log(x) / x

    z = 1.0 / x
    if x < 8.0:                 # y > exp(-32) = 1.2664165549e-14
        x1 = z * _polevl(z, P1, 8) / _p1evl(z, Q1, 8)
    else:
        x1 = z * _polevl(z, P2, 8) / _p1evl(z, Q2, 8)

    x = x0 - x1
    if sign_flag != 0:
        x = -x

    return x
557120
def hessian_local_log_likelihood(self, x):
        
        # Observation likelihoods
        lmbda = np.exp(np.dot(x, self.C.T) + np.dot(self.inputs, self.D.T))
        return np.einsum('tn, ni, nj ->tij', -lmbda, self.C, self.C)
367770
def derivatives(self, x, y, theta_E, s, q):
        
        if q >= 1:
            q = 0.999999
        psi = self._psi(x, y, q, s)
        f_x = theta_E / np.sqrt(1. - q ** 2) * np.arctan(np.sqrt(1. - q ** 2) * x / (psi+s))
        f_y = theta_E / np.sqrt(1. - q ** 2) * np.arctanh(np.sqrt(1. - q ** 2) * y / (psi + q**2*s))
        return f_x, f_y
991485
def log_interp1d(self, xx, yy, kind='linear'):
        
        logx = np.log10(xx)
        logy = np.log10(yy)
        lin_interp = interp1d(logx, logy, kind=kind)
        log_interp = lambda zz: np.power(10.0, lin_interp(np.log10(zz)))
        return log_interp
966653
def ricker(f, length, dt):
    
    t = np.linspace(-int(length/2), int((length-dt)/2), int(length/dt))
    y = (1. - 2.*(np.pi**2)*(f**2)*(t**2))*np.exp(-(np.pi**2)*(f**2)*(t**2))
    return t, y
358188
def derivative(self, t, n=1):
        """"""
        p = self.bpoints()
        if n == 1:
            return 3*(p[1] - p[0])*(1 - t)**2 + 6*(p[2] - p[1])*(1 - t)*t + 3*(
                p[3] - p[2])*t**2
        elif n == 2:
            return 6*(
                (1 - t)*(p[2] - 2*p[1] + p[0]) + t*(p[3] - 2*p[2] + p[1]))
        elif n == 3:
            return 6*(p[3] - 3*(p[2] - p[1]) - p[0])
        elif n > 3:
            return 0
        else:
            raise ValueError("n should be a positive integer.")
1259138
def derivatives(self, x, y, coeffs, beta, center_x=0, center_y=0):
        
        shapelets = self._createShapelet(coeffs)
        r, phi = param_util.cart2polar(x, y, center=np.array([center_x, center_y]))
        alpha1_shapelets, alpha2_shapelets = self._alphaShapelets(shapelets, beta)
        f_x = self._shapeletOutput(r, phi, beta, alpha1_shapelets)
        f_y = self._shapeletOutput(r, phi, beta, alpha2_shapelets)
        return f_x, f_y