Performs a Laplace approximation to the posterior

        Parameters
        ----------
        obj_type : method
            Whether a likelihood or a posterior

        Returns
        ----------
        None (plots posterior)
420853
def _laplace_fit(self,obj_type):
        

        # Get Mode and Inverse Hessian information
        y = self.fit(method='PML',printer=False)

        if y.ihessian is None:
            raise Exception("No Hessian information - Laplace approximation cannot be performed")
        else:

            self.latent_variables.estimation_method = 'Laplace'

            theta, Y, scores, states, states_var, X_names = self._categorize_model_output(self.latent_variables.get_z_values())

            # Change this in future
            try:
                latent_variables_store = self.latent_variables.copy()
            except:
                latent_variables_store = self.latent_variables

            return LaplaceResults(data_name=self.data_name,X_names=X_names,model_name=self.model_name,
                model_type=self.model_type, latent_variables=latent_variables_store,data=Y,index=self.index,
                multivariate_model=self.multivariate_model,objective_object=obj_type, 
                method='Laplace',ihessian=y.ihessian,signal=theta,scores=scores,
                z_hide=self._z_hide,max_lag=self.max_lag,states=states,states_var=states_var)
383666
def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None):
        
        assert mean_function is None, "inference with a mean function not implemented"

        # Compute K
        K = kern.K(X)

        #Find mode
        if self.bad_fhat or self.first_run:
            Ki_f_init = np.zeros_like(Y)
            self.first_run = False
        else:
            Ki_f_init = self._previous_Ki_fhat

        Ki_f_init = np.zeros_like(Y)# FIXME: take this out

        f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata)

        #Compute hessian and other variables at mode
        log_marginal, woodbury_inv, dL_dK, dL_dthetaL = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, kern, Y_metadata)

        self._previous_Ki_fhat = Ki_fhat.copy()
        return Posterior(woodbury_vector=Ki_fhat, woodbury_inv=woodbury_inv, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL}
347963
def _use_gl(objs):
    ''' 

    '''
    from ..models.plots import Plot
    return _any(objs, lambda obj: isinstance(obj, Plot) and obj.output_backend == "webgl")
867232
def _assign_posterior(self):
        

        prior_centers = self.get_centers(self.local_prior)
        posterior_centers = self.get_centers(self.local_posterior_)
        posterior_widths = self.get_widths(self.local_posterior_)
        # linear assignment on centers
        cost = distance.cdist(prior_centers, posterior_centers, 'euclidean')
        _, col_ind = linear_sum_assignment(cost)
        # reorder centers/widths based on cost assignment
        self.set_centers(self.local_posterior_, posterior_centers[col_ind])
        self.set_widths(self.local_posterior_, posterior_widths[col_ind])
        return self
260661
def posterior(self, x, Sigma=None):
        
        n = x.shape[0]
        Sigma = np.eye(self.dim) if Sigma is None else Sigma
        Siginv = inv(Sigma)
        Qpost = inv(self.cov) + n * Siginv
        Sigpost = inv(Qpost)
        mupost = (np.matmul(Siginv, self.mean) + 
                  np.matmu(Siginv, np.sum(x, axis=0)))
        return MvNormal(loc=mupost, cov=Sigpost)
986940
def _parameter_constraints(self, theta_E, gamma, q, phi_G, s_scale):
        
        if theta_E < 0:
            theta_E = 0
        if s_scale < 0.00000001:
            s_scale = 0.00000001
        if gamma < 1.2:
            gamma = 1.2
            theta_E = 0
        if gamma > 2.9:
            gamma = 2.9
            theta_E = 0
        if q < 0.01:
            q = 0.01
            theta_E = 0
        if q > 1:
            q = 1.
            theta_E = 0
        return theta_E, gamma, q, phi_G, s_scale
1083043
def object_type_repr(obj):
    
    if obj is None:
        return 'None'
    elif obj is Ellipsis:
        return 'Ellipsis'
    if obj.__class__.__module__ == '__builtin__':
        name = obj.__class__.__name__
    else:
        name = obj.__class__.__module__ + '.' + obj.__class__.__name__
    return '%s object' % name
1106937
def convenience_calc_log_likelihood(self, params):
        
        shapes, intercepts, betas = self.convenience_split_params(params)

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

        kwargs = {"intercept_params": intercepts,
                  "shape_params": shapes,
                  "ridge": self.ridge,
                  "weights": self.weights}
        log_likelihood = cc.calc_log_likelihood(*args, **kwargs)

        return log_likelihood
976735
def eval_hessian(self, ordered_parameters=[], **parameters):
        
        result = super(HessianObjectiveJacApprox, self).__call__(
            ordered_parameters, **parameters
        )
        num_params = len(self.model.params)
        return [np.broadcast_to(
                    np.zeros_like(comp),
                    (num_params, num_params) + comp.shape
                ) for comp in result]
866801
def _prepare_prior(self, prior, _type):
        '''
        
        '''

        if prior is None and not self.auto_scale:
            prior = self.default_priors.get(term=_type + '_flat')

        if isinstance(prior, Prior):
            prior._auto_scale = False
        else:
            _scale = prior
            prior = self.default_priors.get(term=_type)
            prior.scale = _scale
            if prior.scale is not None:
                prior._auto_scale = False
        return prior
403462
def underlying_likelihood(self, binary_outcomes, modelparams, expparams):
        
        original_mps = modelparams[..., self._orig_mps_slice]
        return self.underlying_model.likelihood(binary_outcomes, original_mps, expparams)
973773
def _optimize_fit(self, obj_type=None, **kwargs):
        

        preopt_search = kwargs.get('preopt_search', True) # If user supplied

        if obj_type == self.neg_loglik:
            method = 'MLE'
        else:
            method = 'PML'

        # Starting values - check to see if model has preoptimize method, if not, simply use default starting values
        if preopt_search is True:
            try:
                phi = self._preoptimize_model(self.latent_variables.get_z_starting_values(), method)
                preoptimized = True
            except:
                phi = self.latent_variables.get_z_starting_values()
                preoptimized = False
        else:
            preoptimized = False
            phi = self.latent_variables.get_z_starting_values()

        phi = kwargs.get('start',phi).copy() # If user supplied

        # Optimize using L-BFGS-B
        p = optimize.minimize(obj_type, phi, method='L-BFGS-B', options={'gtol': 1e-8})
        if preoptimized is True:
            p2 = optimize.minimize(obj_type, self.latent_variables.get_z_starting_values(), method='L-BFGS-B', 
                options={'gtol': 1e-8})
            if self.neg_loglik(p2.x) < self.neg_loglik(p.x):
                p = p2

        theta, Y, scores, states, states_var, X_names = self._categorize_model_output(p.x)

        # Check that matrix is non-singular; act accordingly
        try:
            ihessian = np.linalg.inv(nd.Hessian(obj_type)(p.x))
            ses = np.power(np.abs(np.diag(ihessian)),0.5)
            self.latent_variables.set_z_values(p.x,method,ses,None)

        except:
            ihessian = None
            ses = None
            self.latent_variables.set_z_values(p.x,method,None,None)

        self.latent_variables.estimation_method = method

        # Change this in future
        try:
            latent_variables_store = self.latent_variables.copy()
        except:
            latent_variables_store = self.latent_variables

        return MLEResults(data_name=self.data_name,X_names=X_names,model_name=self.model_name,
                model_type=self.model_type, latent_variables=latent_variables_store,results=p,data=Y, index=self.index,
                multivariate_model=self.multivariate_model,objective_object=obj_type, 
                method=method,ihessian=ihessian,signal=theta,scores=scores,
                z_hide=self._z_hide,max_lag=self.max_lag,states=states,states_var=states_var)
919585
def as_ompenv(cls, obj):
        """"""
        if isinstance(obj, cls): return obj
        if obj is None: return cls()
        return cls(**obj)
1247816
def extend_dirichlet(p):
    
    if len(np.shape(p)) > 1:
        return np.hstack((p, np.atleast_2d(1. - np.sum(p))))
    else:
        return np.hstack((p, 1. - np.sum(p)))
955951
def map_parameter_estimate(self):
        
        res = self.pst.res
        assert res is not None
        # build the prior expectation parameter vector
        prior_expt = self.pst.parameter_data.loc[:,["parval1"]].copy()
        islog = self.pst.parameter_data.partrans == "log"
        prior_expt.loc[islog] = prior_expt.loc[islog].apply(np.log10)
        prior_expt = Matrix.from_dataframe(prior_expt)
        prior_expt.col_names = ["prior_expt"]
        # build the residual vector
        res_vec = Matrix.from_dataframe(res.loc[:,["residual"]])

        # form the terms of Schur's complement
        b = self.parcov * self.jco.T
        c = ((self.jco * self.parcov * self.jco.T) + self.obscov).inv
        bc = Matrix((b * c).x, row_names=b.row_names, col_names=c.col_names)

        # calc posterior expectation
        upgrade = bc * res_vec
        upgrade.col_names = ["prior_expt"]
        post_expt = prior_expt + upgrade

        # post processing - back log transform
        post_expt = pd.DataFrame(data=post_expt.x,index=post_expt.row_names,
                                 columns=["post_expt"])
        post_expt.loc[:,"prior_expt"] = prior_expt.x.flatten()
        post_expt.loc[islog,:] = 10.0**post_expt.loc[islog,:]
        # unc_sum = self.get_parameter_summary()
        # post_expt.loc[:,"standard_deviation"] = unc_sum.post_var.apply(np.sqrt)
        post_expt.sort_index(inplace=True)
        return post_expt
988236
def _laplace_approximation_newton(self, tol=1e-6, stepsz=0.9, verbose=False):
        
        from pylds.util import solve_symm_block_tridiag, scipy_solve_symm_block_tridiag
        scale = self.T * self.D_emission

        def newton_step(x, stepsz):
            assert 0 <= stepsz <= 1
            g = self.gradient_log_joint(x)
            H_diag, H_upper_diag = self.sparse_hessian_log_joint(x)
            Hinv_g = -scipy_solve_symm_block_tridiag(-H_diag / scale,
                                               -H_upper_diag / scale,
                                               g / scale)
            return x - stepsz * Hinv_g

        if verbose:
            print("Fitting Laplace approximation")

        itr = [0]
        def cbk(x):
            print("Iteration: ", itr[0],
                  "\tObjective: ", (self.log_joint(x) / scale).round(4),
                  "\tAvg Grad: ", (self.gradient_log_joint(x).mean() / scale).round(4))
            itr[0] += 1

        # Solve for optimal x with Newton's method
        x = self.gaussian_states
        dx = np.inf
        while dx >= tol:
            xnew = newton_step(x, stepsz)
            dx = np.mean(abs(xnew - x))
            x = xnew

            if verbose:
                cbk(x)

        assert np.all(np.isfinite(x))
        if verbose:
            print("Done")

        return x
875960
def _bool_method_SERIES(cls, op, special):
    
    op_name = _get_op_name(op, special)

    def na_op(x, y):
        try:
            result = op(x, y)
        except TypeError:
            assert not isinstance(y, (list, ABCSeries, ABCIndexClass))
            if isinstance(y, np.ndarray):
                # bool-bool dtype operations should be OK, should not get here
                assert not (is_bool_dtype(x) and is_bool_dtype(y))
                x = ensure_object(x)
                y = ensure_object(y)
                result = libops.vec_binop(x, y, op)
            else:
                # let null fall thru
                assert lib.is_scalar(y)
                if not isna(y):
                    y = bool(y)
                try:
                    result = libops.scalar_binop(x, y, op)
                except (TypeError, ValueError, AttributeError,
                        OverflowError, NotImplementedError):
                    raise TypeError("cannot compare a dtyped [{dtype}] array "
                                    "with a scalar of type [{typ}]"
                                    .format(dtype=x.dtype,
                                            typ=type(y).__name__))

        return result

    fill_int = lambda x: x.fillna(0)
    fill_bool = lambda x: x.fillna(False).astype(bool)

    def wrapper(self, other):
        is_self_int_dtype = is_integer_dtype(self.dtype)

        self, other = _align_method_SERIES(self, other, align_asobject=True)
        res_name = get_op_result_name(self, other)

        if isinstance(other, ABCDataFrame):
            # Defer to DataFrame implementation; fail early
            return NotImplemented

        elif isinstance(other, (ABCSeries, ABCIndexClass)):
            is_other_int_dtype = is_integer_dtype(other.dtype)
            other = fill_int(other) if is_other_int_dtype else fill_bool(other)

            ovalues = other.values
            finalizer = lambda x: x

        else:
            # scalars, list, tuple, np.array
            is_other_int_dtype = is_integer_dtype(np.asarray(other))
            if is_list_like(other) and not isinstance(other, np.ndarray):
                # TODO: Can we do this before the is_integer_dtype check?
                # could the is_integer_dtype check be checking the wrong
                # thing?  e.g. other = [[0, 1], [2, 3], [4, 5]]?
                other = construct_1d_object_array_from_listlike(other)

            ovalues = other
            finalizer = lambda x: x.__finalize__(self)

        # For int vs int `^`, `|`, `&` are bitwise operators and return
        #   integer dtypes.  Otherwise these are boolean ops
        filler = (fill_int if is_self_int_dtype and is_other_int_dtype
                  else fill_bool)
        res_values = na_op(self.values, ovalues)
        unfilled = self._constructor(res_values,
                                     index=self.index, name=res_name)
        filled = filler(unfilled)
        return finalizer(filled)

    wrapper.__name__ = op_name
    return wrapper