Calculates and returns the negative of the log-likelihood and the
        negative of the gradient. This function is used as the objective
        function in scipy.optimize.minimize.
964807
def calc_neg_log_likelihood_and_neg_gradient(self, params):
        
        neg_log_likelihood = -1 * self.convenience_calc_log_likelihood(params)
        neg_gradient = -1 * self.convenience_calc_gradient(params)

        if self.constrained_pos is not None:
            neg_gradient[self.constrained_pos] = 0

        return neg_log_likelihood, neg_gradient
990437
def _extract_log_likelihood(self):
        
        # This next line is brittle and may not work forever, but is a secret
        # way to access the model from the trace.
        model = self.trace._straces[0].model  # pylint: disable=protected-access
        if len(model.observed_RVs) != 1:
            return None, None
        else:
            if self.dims is not None:
                coord_name = self.dims.get(model.observed_RVs[0].name)
            else:
                coord_name = None

        cached = [(var, var.logp_elemwise) for var in model.observed_RVs]

        def log_likelihood_vals_point(point):
            """Compute log likelihood for each observed point."""
            log_like_vals = []
            for var, log_like in cached:
                log_like_val = log_like(point)
                if var.missing_values:
                    log_like_val = log_like_val[~var.observations.mask]
                log_like_vals.append(log_like_val.ravel())
            return np.concatenate(log_like_vals)

        chain_likelihoods = []
        for chain in self.trace.chains:
            log_like = (log_likelihood_vals_point(point) for point in self.trace.points([chain]))
            chain_likelihoods.append(np.stack(log_like))
        return np.stack(chain_likelihoods), coord_name
247691
def neg_loglik(self, beta):
             
        Z = np.zeros(2)
        Z[0] = 1          
        states = np.zeros([self.state_no, self.data.shape[0]])
        states[0,:] = beta[self.z_no:self.z_no+self.data.shape[0]] 
        states[1,:] = beta[self.z_no+self.data.shape[0]:] 
        parm = np.array([self.latent_variables.z_list[k].prior.transform(beta[k]) for k in range(self.z_no)]) # transformed distribution parameters
        scale, shape, skewness = self._get_scale_and_shape(parm)
        return self.state_likelihood(beta, states) + self.family.neg_loglikelihood(self.data, self.link(np.dot(Z, states)), scale, shape, skewness)
771625
def logp_partial_gradient(self, variable, calculation_set=None):
        
        if self.verbose > 0:
            print_('\t' + self.__name__ + ': logp_partial_gradient accessed.')

        if not (datatypes.is_continuous(variable)
                and datatypes.is_continuous(self)):
            return zeros(shape(variable.value))

        # loop through all the parameters and add up all the gradients of log p
        # with respect to the approrpiate variable
        gradient = builtins.sum(
            [child.logp_partial_gradient(self,
                                         calculation_set) for child in self.children])

        totalGradient = 0
        for parameter, value in six.iteritems(self.parents):
            if value is variable:

                totalGradient += self.apply_jacobian(
                    parameter, variable, gradient)

        return np.reshape(totalGradient, shape(variable.value))
774664
def likelihood(self):
        '''
        
        '''
        try:
            likelihood = self.__success / (self.__success + self.__failure)
        except ZeroDivisionError:
            likelihood = 0.0
        return likelihood
86235
def minimise(objective_function, p0, **kwargs):
    

    p0 = np.array(p0)

    op_kwargs = kwargs.copy()

    # Which optimisation algorithm?
    available_methods = ("BFGS", "Nelder-Mead", "Powell", "CG", "TNC")
    method = op_kwargs.pop("method", available_methods[0])
    logger.debug("Optimising with {0} algorithm".format(method))

    op_kwargs.update({ "callback": None, "full_output": True, "retall": False })

    def _drop_forbidden_keywords(op_kwds, allowed):
        for keyword in set(op_kwds).difference(allowed):
            logger.debug("Dropping optimisation keyword {0}: {1}".format(
                keyword, op_kwds.pop(keyword)))

    t_init = time()
    if method == "Nelder-Mead":
        # Need to force some arguments.
        # Optional arguments:
        # xtol, ftol, maxiter, maxfun
        _drop_forbidden_keywords(op_kwargs, ("xtol", "ftol", "maxiter", 
            "maxfun", "full_output", "retall", "callback", "disp"))

        # Set defaults.
        op_kwargs.setdefault("disp", False)

        x_opt, f_opt, num_iter, num_funcalls, warnflag \
            = op.fmin(objective_function, p0, **op_kwargs)

        logger.debug("Number of iterations: {0}, number of function calls: {1}"\
            .format(num_iter, num_funcalls))

        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.warn("Optimisation failed: {0}".format([
                "Maximum number of function evaluations.",
                "Maximum number of iterations."][warnflag - 1]))

    elif method == "Powell":
        # Optional arguments: xtol, ftol, maxiter, maxfun
        _drop_forbidden_keywords(op_kwargs, ("xtol", "ftol", "maxiter", 
            "maxfun", "retall", "callback", "full_output", "disp"))
        
        # Set defaults.
        op_kwargs.setdefault("disp", False)

        x_opt, f_opt, direc, num_iter, num_funcalls, warnflag \
            = op.fmin_powell(objective_function, p0, **op_kwargs)

        logger.debug("Number of iterations: {0}, number of function calls: {1}"\
            .format(num_iter, num_funcalls))

        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.warn("Optimisation failed: {0}".format([
                "Maximum number of function evaluations.",
                "Maximum number of iterations."][warnflag - 1]))

    elif method == "CG":
        # Optional arguments: gtol, norm, epsilon, maxiter.
        _drop_forbidden_keywords(op_kwargs, ("gtol", "norm", "epsilon", 
            "maxiter", "retall", "callback", "full_output"))

        x_opt, f_opt, num_funcalls, num_gradcalls, warnflag \
            = op.fmin_cg(objective_function, p0, **op_kwargs)

        logger.debug("Number of function calls: {0}, gradient calls: {1}"\
            .format(num_funcalls, num_gradcalls))

        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.debug("Optimisation failed: {0}".format([
                "Maximum number of iterations exceeded",
                "Gradient and/or function calls were not changing."
                ][warnflag - 1]))

    elif method == "BFGS":
        # Since we have at least some boundaries, this will actually call
        # the L-BFGS-B algorithm.
        # Optional arguments: m, factr, pgtol, epsilon, maxfun, maxiter

        _drop_forbidden_keywords(op_kwargs, ("m", "factr", "pgtol", 
            "approx_grad", "epsilon", "disp", "iprint", "maxfun", "maxiter"))

        # Default/required:
        op_kwargs.setdefault("factr", 10.0)
        op_kwargs["approx_grad"] = True
        
        # Because the parameters will vary by orders of magnitude, here we
        # scale everything to the initial value so that the epsilon keyword
        # makes some sense.

        scale = p0.copy()
        
        def scaled_objective_function(theta):
            return objective_function(theta.copy() * scale)

        x_opt, f_opt, info_dict = op.fmin_l_bfgs_b(scaled_objective_function,
            np.ones(p0.size, dtype=float), **op_kwargs)

        # Rescale.
        x_opt *= scale

        g_opt = info_dict["grad"]
        num_iter = info_dict["nit"]
        num_funcalls = info_dict["funcalls"]
        warnflag = info_dict["warnflag"]

        logger.debug("Number of iterations: {0}, function calls: {1}".format(
            num_iter, num_funcalls))
        if warnflag == 0:
            logger.debug("Optimisation converged after {0:.1f} seconds.".format(
                time() - t_init))
        else:
            logger.warn("Optimisation failed: {0}".format([
                "Too many function evaluations or iterations!",
                "{0}".format(info_dict.get("task", None))][warnflag - 1]))
    
    elif method == "TNC":
        # Optional arguments:   bounds, epsilon, scale, offset, messages, 
        #                       maxCGit, maxfun, eta, stepmx, accuracy, 
        #                       fmin, ftol, xtol, pgtol, rescale

        _drop_forbidden_keywords(op_kwargs, ("approx_grad", "bounds", "epsilon",
            "scale", "offset", "messages", "maxCGit", "maxfun", "eta", "stepmx",
            "accuracy", "fmin", "ftol", "xtol", "pgtol", "rescale", "disp"))

        # Required:
        op_kwargs["approx_grad"] = True

        x_opt, num_funcalls, rc = op.fmin_tnc(objective_function, p0, 
            **op_kwargs)
        
        rcstring = {
            -1: "Infeasible (lower bound > upper bound)",
            0: "Local minimum reached (|pg| ~= 0)",
            1: "Converged (|f_n-f_(n-1)| ~= 0)",
            2: "Converged (|x_n-x_(n-1)| ~= 0)",
            3: "Max. number of function evaluations reached",
            4: "Linear search failed",
            5: "All lower bounds are equal to the upper bounds",
            6: "Unable to progress",
            7: "User requested end of minimization"
        }[rc]

        logger.debug("Number of function calls: {0}, result: {1}".format(
            num_funcalls, rcstring))
        if rc in (1, 2):
            logger.debug("Optimisation converged after {0:.1f} seconds: {1}"\
                .format(time() - t_init, rcstring))
        else:
            logger.warn("Optimisation failed: {0}".format(rcstring))

    else:
        raise ValueError("optimisation algorithm {0} is not available "\
            "(available methods are {1})".format(
                method, ", ".join(available_methods)))


    return x_opt
329466
def convenience_calc_log_likelihood(self, params):
        
        shapes, intercepts, betas = self.convenience_split_params(params)

        args = [betas,
                self.design_3d,
                self.alt_id_vector,
                self.rows_to_obs,
                self.rows_to_alts,
                self.rows_to_mixers,
                self.choice_vector,
                self.utility_transform]

        kwargs = {"ridge": self.ridge, "weights": self.weights}
        log_likelihood = general_log_likelihood(*args, **kwargs)

        return log_likelihood
197406
def fit(self, x0=None, distribution='lognormal', n=None, **kwargs):
        '''
        '''
        dist = {'lognormal': PSDLognormal, 
                'GGS': PSDGatesGaudinSchuhman, 
                'RR': PSDRosinRammler}[distribution]
        
        if distribution == 'lognormal':
            if x0 is None:
                d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                s = 0.4
                x0 = [d_characteristic, s]
        elif distribution == 'GGS':
            if x0 is None:
                d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                m = 1.5
                x0 = [d_characteristic, m]
        elif distribution == 'RR':
            if x0 is None:
                x0 = [5E-6, 1e-2]
        from scipy.optimize import minimize
        return minimize(self._fit_obj_function, x0, args=(dist, n), **kwargs)
552008
def _fit(
        self,
        minimizing_function_args,
        iterative_fitting,
        initial_params,
        params_size,
        disp,
        tol=1e-6,
        fit_method="Nelder-Mead",
        maxiter=2000,
        **kwargs
    ):
        """"""
        ll = []
        sols = []

        if iterative_fitting <= 0:
            raise ValueError("iterative_fitting parameter should be greater than 0 as of lifetimes v0.2.1")

        if iterative_fitting > 1 and initial_params is not None:
            raise ValueError(
                "iterative_fitting and initial_params should not be both set, as no improvement could be made."
            )

        # set options for minimize, if specified in kwargs will be overwritten
        minimize_options = {}
        minimize_options["disp"] = disp
        minimize_options["maxiter"] = maxiter
        minimize_options.update(kwargs)

        total_count = 0
        while total_count < iterative_fitting:
            current_init_params = (
                np.random.normal(1.0, scale=0.05, size=params_size) if initial_params is None else initial_params
            )
            if minimize_options["disp"]:
                print("Optimize function with {}".format(fit_method))

            output = minimize(
                self._negative_log_likelihood,
                method=fit_method,
                tol=tol,
                x0=current_init_params,
                args=minimizing_function_args,
                options=minimize_options,
            )
            sols.append(output.x)
            ll.append(output.fun)

            total_count += 1
        argmin_ll, min_ll = min(enumerate(ll), key=lambda x: x[1])
        minimizing_params = sols[argmin_ll]
        return minimizing_params, min_ll
612487
def brent_optimise(node1, node2, min_brlen=0.001, max_brlen=10, verbose=False):
    
    from scipy.optimize import minimize_scalar
    wrapper = BranchLengthOptimiser(node1, node2, (min_brlen + max_brlen) / 2.)
    n = minimize_scalar(lambda x: -wrapper(x)[0], method='brent', bracket=(min_brlen, max_brlen))['x']
    if verbose:
        logger.info(wrapper)
    if n < min_brlen:
        n = min_brlen
        wrapper(n)
    return n, -1 / wrapper.get_d2lnl(n)
192088
def p_expr_minus_expr(p):
    
    p[0] = make_binary(p.lineno(2), 'MINUS', p[1], p[3], lambda x, y: x - y)
1187803
def hessian_local_log_likelihood(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        hfun = hessian(self.local_log_likelihood)
        H_diag = np.zeros((T, D, D))
        for t in range(T):
            H_diag[t] = hfun(x[t], self.data[t], self.inputs[t])
        return H_diag
262307
def get_err(snr):
    

    quad_terms = np.array(
            [3.11e-3, 1.10e-5, 6.95e-6, 5.05e-6, 4.65e-6, 4.10e-6])
    lin_terms = np.array(
            [-0.869, -2.07e-3, -1.40e-3, -1.03e-3, -1.13e-3, -7.29e-4])
    consts = np.array([104, 0.200, 0.117, 0.114, 0.156, 0.0624])
    err = quad_terms[:,None] * snr**2 + lin_terms[:,None] * snr + consts[:,None]

    # find the minimum of the quadratic function
    min_snr = -lin_terms / (2*quad_terms)
    min_err = quad_terms * min_snr**2 + lin_terms * min_snr + consts
    mask = (snr[:,None] > min_snr).T
    for i in range(0,len(min_err)):
        err[i][mask[i]] = min_err[i]

    return err
1329345
def _build_likelihood(self):
        

        # FITC approximation to the log marginal likelihood is
        # log ( normal( y | mean, K_fitc ) )
        # where K_fitc = Qff + diag( \nu )
        # where Qff = Kfu Kuu^{-1} Kuf
        # with \nu_i = Kff_{i,i} - Qff_{i,i} + \sigma^2

        # We need to compute the Mahalanobis term -0.5* err^T K_fitc^{-1} err
        # (summed over functions).

        # We need to deal with the matrix inverse term.
        # K_fitc^{-1} = ( Qff + \diag( \nu ) )^{-1}
        #            = ( V^T V + \diag( \nu ) )^{-1}
        # Applying the Woodbury identity we obtain
        #            = \diag( \nu^{-1} ) - \diag( \nu^{-1} ) V^T ( I + V \diag( \nu^{-1} ) V^T )^{-1) V \diag(\nu^{-1} )
        # Let \beta =  \diag( \nu^{-1} ) err
        # and let \alpha = V \beta
        # then Mahalanobis term = -0.5* ( \beta^T err - \alpha^T Solve( I + V \diag( \nu^{-1} ) V^T, alpha ) )

        err, nu, Luu, L, alpha, beta, gamma = self._build_common_terms()

        mahalanobisTerm = -0.5 * tf.reduce_sum(tf.square(err) / tf.expand_dims(nu, 1)) \
                          + 0.5 * tf.reduce_sum(tf.square(gamma))

        # We need to compute the log normalizing term -N/2 \log 2 pi - 0.5 \log \det( K_fitc )

        # We need to deal with the log determinant term.
        # \log \det( K_fitc ) = \log \det( Qff + \diag( \nu ) )
        #                    = \log \det( V^T V + \diag( \nu ) )
        # Applying the determinant lemma we obtain
        #                    = \log [ \det \diag( \nu ) \det( I + V \diag( \nu^{-1} ) V^T ) ]
        #                    = \log [ \det \diag( \nu ) ] + \log [ \det( I + V \diag( \nu^{-1} ) V^T ) ]

        constantTerm = -0.5 * self.num_data * tf.log(tf.constant(2. * np.pi, settings.float_type))
        logDeterminantTerm = -0.5 * tf.reduce_sum(tf.log(nu)) - tf.reduce_sum(tf.log(tf.matrix_diag_part(L)))
        logNormalizingTerm = constantTerm + logDeterminantTerm

        return mahalanobisTerm + logNormalizingTerm * self.num_latent
37917
def log_loss(preds, labels):
    """"""
    log_likelihood = np.sum(labels * np.log(preds)) / len(preds)
    return -log_likelihood
1208996
function(q, lowerTail, logp) {
            var ind, ret;

            logp = logp === true;
            lowerTail = lowerTail !== false;

            if (utils.hasNaN(q)) { return NaN; }
            if (q < xs[0]) {
               ret = lowerTail ? 0 : 1;
            } else if (q >= xs[xs.length - 1]) {
               ret = lowerTail ? 1 : 0;
            } else {
               ind = binSearch(xs, function(y) { return y <= q; });
               ret = lowerTail ? cumsLeft[ind] : cumsRight[ind];
            }
            ret = Math.min(ret, 1);
            return logp ? Math.log(ret) : ret;
         }
474259
def (z):
    
    if z < -6:
        return 1, -1.0e12, -1
    if z > 6:
        return 0, 0, 1
    else:
        logphi = -0.5 * (z * z + l2p)
        logPhi = np.log(.5 * special.erfc(-z / sq2))
        e = np.exp(logphi - logPhi)
    return e, logPhi, 0