One dimensional power law with an exponential cutoff model function
180485
def eval(e, amplitude, e_0, alpha, e_cutoff, beta):
        

        xx = e / e_0
        return amplitude * xx ** (-alpha) * np.exp(-(e / e_cutoff) ** beta)
542295
def _expectation(p, kern, none1, none2, none3, nghp=None):
    
    return kern.Kdiag(p.mu)
614325
function ecp_dbl(R, N, a, p, l) {
	var x = N[0], y = N[1], z = N[2];
	var xr = R[0], yr = R[1], zr = R[2];
	if (mpn_zero_p(z, l)) {
		mpn_zero(xr, l), mpn_zero(yr, l), mpn_zero(zr, l);
		return R;
	}
	var t0 = _ecp_t0, t1 = _ecp_t1, t2 = _ecp_t2, t3 = _ecp_t3;
	fp_add(t0, t0, fp_dbl(t1, fp_sqr(t0, x, p, l), p, l), p, l);
	if (!mpn_zero_p(a, l)) {
		fp_add(t0, t0, fp_mul(t1, a, fp_sqr(t2, fp_sqr(t1, z, p, l), p, l), p, l), p, l);
	}
	fp_dbl(t1, fp_sqr(t1, y, p, l), p, l);
	fp_dbl(t2, fp_mul(t2, x, t1, p, l), p, l);
	fp_dbl(t3, fp_sqr(t3, t1, p, l), p, l);
	fp_sub(xr, fp_sqr(xr, t0, p, l), fp_dbl(t1, t2, p, l), p, l);
	fp_sub(yr, fp_mul(yr, t0, fp_sub(t1, t2, xr, p, l), p, l), t3, p, l);
	fp_dbl(zr, fp_mul(zr, y, z, p, l), p, l);
	return R;
}
1059980
def hessian(self, x, y, kappa_ext, ra_0=0, dec_0=0):
        
        gamma1 = 0
        gamma2 = 0
        kappa = kappa_ext
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy
884310
def _twofilter_smoothing_ON(self, t, ti, info, phi, lwinfo, return_ess,
                               modif_forward, modif_info):
        
        if modif_info is not None:
            lwinfo += modif_info
        Winfo = rs.exp_and_normalise(lwinfo)
        I = rs.multinomial(Winfo)
        if modif_forward is not None:
            lw = self.wgt[t].lw + modif_forward
            W = rs.exp_and_normalise(lw)
        else:
            W = self.wgt[t].W
        J = rs.multinomial(W)
        log_omega = self.model.logpt(t + 1, self.X[t][J], info.hist.X[ti][I])
        if modif_forward is not None:
            log_omega -= modif_forward[J]
        if modif_info is not None:
            log_omega -= modif_info[I]
        Om = rs.exp_and_normalise(log_omega)
        est = np.average(phi(self.X[t][J], info.hist.X[ti][I]), axis=0,
                         weights=Om)
        if return_ess:
            return (est, 1. / np.sum(Om**2))
        else:
            return est
394987
def _M(self, k, t, tips=None, gaps=None):
        """"""
        if self._distributionmodel:
            return self.model.M(k, t, tips, gaps)
        else:
            return self.model.M(t, tips, gaps)
762343
def estimate(self, maxiter=250, convergence=1e-7):

        

        self.loglik = np.zeros(maxiter)

        iter = 0

        while iter < maxiter:

            self.loglik[iter] = self.E_step()
            if np.isnan(self.loglik[iter]):
                    print("undefined log-likelihood")
                    break
            self.M_step()

            if self.loglik[iter] - self.loglik[iter - 1] < 0 and iter > 0:
                    print("log-likelihood decreased by %f at iteration %d"
                          % (self.loglik[iter] - self.loglik[iter - 1],
                             iter))
            elif self.loglik[iter] - self.loglik[iter - 1] < convergence \
                    and iter > 0:
                    print("convergence at iteration %d, loglik = %f" %
                           (iter, self.loglik[iter]))
                    self.loglik = self.loglik[self.loglik < 0]
                    break

            iter += 1
121592
def parametric_mean_function(max_iters=100, optimize=True, plot=True):
    
    #create  simple mean function
    mf = GPy.core.Mapping(1,1)
    mf.f = np.sin

    X = np.linspace(0,10,50).reshape(-1,1)
    Y = np.sin(X) + 0.5*np.cos(3*X) + 0.1*np.random.randn(*X.shape) + 3*X

    mf = GPy.mappings.Linear(1,1)

    k =GPy.kern.RBF(1)
    lik = GPy.likelihoods.Gaussian()
    m = GPy.core.GP(X, Y, kernel=k, likelihood=lik, mean_function=mf)
    if optimize:
        m.optimize(max_iters=max_iters)
    if plot:
        m.plot()
    return m
965549
def toy_poisson_rbf_1d_laplace(optimize=True, plot=True):
    """"""
    optimizer='scg'
    x_len = 100
    X = np.linspace(0, 10, x_len)[:, None]
    f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.RBF(1).K(X))
    Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:,None]

    kern = GPy.kern.RBF(1)
    poisson_lik = GPy.likelihoods.Poisson()
    laplace_inf = GPy.inference.latent_function_inference.Laplace()

    # create simple GP Model
    m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf)

    if optimize:
        m.optimize(optimizer)
    if plot:
        m.plot()
        # plot the real underlying rate function
        pb.plot(X, np.exp(f_true), '--k', linewidth=2)

    return m
38821
def power_law_anisotropy(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        # first term
        theta_E = kwargs_profile['theta_E']
        gamma = kwargs_profile['gamma']
        r_ani = kwargs_anisotropy['r_ani']
        a = 0.551 * kwargs_light['r_eff']
        rho0_r0_gamma = self._rho0_r0_gamma(theta_E, gamma)
        prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
        prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
        hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
        hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
        fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
        sigma2_dim_less = prefac1 * prefac2 * fac
        return sigma2_dim_less * (self.cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000)**2
711236
def nhill_i(self,x,threshold=0.1,power=2):
        
        x_pow = np.power(x,power)
        threshold_pow = np.power(threshold,power)
        return threshold_pow / (x_pow + threshold_pow) * (1 - x_pow)
66008
def hessian_local_log_likelihood(self, x):
        
        C, D, u, y = self.C, self.D, self.inputs, self.data
        psi = x.dot(C.T) + u.dot(D.T)
        p = 1. / (1 + np.exp(-psi))
        dp_dpsi = p * (1 - p)
        return np.einsum('tn, ni, nj ->tij', -dp_dpsi, self.C, self.C)
728607
def sse(mean, estimator):
    
    return np.sum((np.asarray(estimator) - np.asarray(mean)) ** 2, axis=0)
218408
def kld(d1, d2):
    
    d1, d2 = flatten(d1), flatten(d2)
    return entropy(d1, d2, 2.0)
707982
def F_t(X, Y, S, M_E, E, m0, rho):
    '''  '''
    r = X.shape[1]
    out1 = (((np.dot(np.dot(X, S), Y.T) - M_E) * E)**2).sum() / 2
    out2 = rho * G(Y, m0, r)
    out3 = rho * G(X, m0, r)

    return out1 + out2 + out3
59082
def rv_model(ts, K, e, omega, chi, P):
    """"""

    assert ts.ndim == 1, 'ts must be one-dimensional'

    return kp.rv_model(ts, K, e, omega, chi, P)
331297
def _expectation(p, kern, none1, none2, none3, nghp=None):
    
    with params_as_tensors_for(kern):
        # use only active dimensions
        Xmu, _ = kern._slice(p.mu, None)
        Xcov = kern._slice_cov(p.cov)

        return tf.reduce_sum(kern.variance * (tf.matrix_diag_part(Xcov) + tf.square(Xmu)), 1)