def calc_bhhh_hessian_approximation_mixed_logit(params,
                                                design_3d,
                                                alt_IDs,
                                                rows_to_obs,
                                                rows_to_alts,
                                                rows_to_mixers,
                                                choice_vector,
                                                utility_transform,
                                                ridge=None,
                                                weights=None):
    """
    Parameters
    ----------
    params : 1D ndarray.
        All elements should by ints, floats, or longs. Should have 1 element
        for each utility coefficient being estimated (i.e. num_features  +
        num_coefs_being_mixed).
    design_3d : 3D ndarray.
        All elements should be ints, floats, or longs. Should have one row per
        observation per available alternative. The second axis should have as
        many elements as there are draws from the mixing distributions of the
        coefficients. The last axis should have one element per index
        coefficient being estimated.
    alt_IDs : 1D ndarray.
        All elements should be ints. There should be one row per obervation per
        available alternative for the given observation. Elements denote the
        alternative corresponding to the given row of the design matrix.
    rows_to_obs : 2D scipy sparse array.
        All elements should be zeros and ones. Should have one row per
        observation per available alternative and one column per observation.
        This matrix maps the rows of the design matrix to the unique
        observations (on the columns).
    rows_to_alts : 2D scipy sparse array.
        All elements should be zeros and ones. Should have one row per
        observation per available alternative and one column per possible
        alternative. This matrix maps the rows of the design matrix to the
        possible alternatives for this dataset.
    rows_to_mixers : 2D scipy sparse array.
        All elements should be zeros and ones. Will map the rows of the design
        matrix to the particular units that the mixing is being performed over.
        Note that in the case of panel data, this matrix will be different from
        `rows_to_obs`.
    choice_vector : 1D ndarray.
        All elements should be either ones or zeros. There should be one row
        per observation per available alternative for the given observation.
        Elements denote the alternative which is chosen by the given
        observation with a 1 and a zero otherwise.
    utility_transform :  callable.
        Should accept a 1D array of systematic utility values, a 1D array of
        alternative IDs, and miscellaneous args and kwargs. Should return a 2D
        array whose elements contain the appropriately transformed systematic
        utility values, based on the current model being evaluated and the
        given draw of the random coefficients. There should be one column for
        each draw of the random coefficients. There should have one row per
        individual per choice situation per available alternative.
    ridge : int, float, long, or None, optional.
        Determines whether or not ridge regression is performed. If a float is
        passed, then that float determines the ridge penalty for the
        optimization. Default = None.
    weights : 1D ndarray or None, optional.
        Allows for the calculation of weighted log-likelihoods. The weights can
        represent various things. In stratified samples, the weights may be
        the proportion of the observations in a given strata for a sample in
        relation to the proportion of observations in that strata in the
        population. In latent class models, the weights may be the probability
        of being a particular class. Default == None.

    Returns
    -------
    bhhh_matrix : 2D ndarray of shape `(design.shape[1], design.shape[1])`.
        The returned array is the BHHH approximation of the Fisher Information
        Matrix. I.e it is the negative of the sum of the outer product of
        each individual's gradient with itself.
    """
    # Calculate the weights for the sample
    if weights is None:
        weights = np.ones(design_3d.shape[0])
    weights_per_obs =\
        np.max(rows_to_mixers.toarray() * weights[:, None], axis=0)
    # Calculate the regular probability array. Note the implicit assumption
    # that params == index coefficients.
    prob_array = general_calc_probabilities(params,
                                            design_3d,
                                            alt_IDs,
                                            rows_to_obs,
                                            rows_to_alts,
                                            utility_transform,
                                            return_long_probs=True)

    # Calculate the simulated probability of correctly predicting each persons
    # sequence of choices. Note that this function implicitly assumes that the
    # mixing unit is the individual
    prob_results = calc_choice_sequence_probs(prob_array,
                                              choice_vector,
                                              rows_to_mixers,
                                              return_type="all")
    # Calculate the sequence probabilities given random draws
    # and calculate the overal simulated probabilities
    sequence_prob_array = prob_results[1]
    simulated_probs = prob_results[0]

    # Convert the various probabilties to long format
    long_sequence_prob_array = rows_to_mixers.dot(sequence_prob_array)
    long_simulated_probs = rows_to_mixers.dot(simulated_probs)
    # Scale sequence probabilites given random draws by simulated probabilities
    scaled_sequence_probs = (long_sequence_prob_array /
                             long_simulated_probs[:, None])
    # Calculate the scaled error. Will have shape == (num_rows, num_draws)
    scaled_error = ((choice_vector[:, None] - prob_array) *
                    scaled_sequence_probs)

    # Calculate the gradient. Note that the lines below assume that we are
    # taking the gradient of an MNL model. Should refactor to make use of the
    # built in gradient function for logit-type models. Should also refactor
    # the gradient function for logit-type models to be able to handle 2D
    # systematic utility arrays. `gradient` will have shape
    # (design_3d.shape[0], design_3d.shape[2])
    gradient = (scaled_error[:, :, None] * design_3d).mean(axis=1)

    gradient_per_obs = rows_to_mixers.T.dot(gradient)

    bhhh_matrix =\
        gradient_per_obs.T.dot(weights_per_obs[:, None] * gradient_per_obs)

    if ridge is not None:
        bhhh_matrix -= 2 * ridge

    # Note the "-1" is because we are approximating the Fisher information
    # matrix which has a negative one in the front of it?
    return -1 * bhhh_matrix
Parameters
    ----------
    params : 1D ndarray.
        All elements should by ints, floats, or longs. Should have 1 element
        for each utility coefficient being estimated (i.e. num_features  +
        num_coefs_being_mixed).
    design_3d : 3D ndarray.
        All elements should be ints, floats, or longs. Should have one row per
        observation per available alternative. The second axis should have as
        many elements as there are draws from the mixing distributions of the
        coefficients. The last axis should have one element per index
        coefficient being estimated.
    alt_IDs : 1D ndarray.
        All elements should be ints. There should be one row per obervation per
        available alternative for the given observation. Elements denote the
        alternative corresponding to the given row of the design matrix.
    rows_to_obs : 2D scipy sparse array.
        All elements should be zeros and ones. Should have one row per
        observation per available alternative and one column per observation.
        This matrix maps the rows of the design matrix to the unique
        observations (on the columns).
    rows_to_alts : 2D scipy sparse array.
        All elements should be zeros and ones. Should have one row per
        observation per available alternative and one column per possible
        alternative. This matrix maps the rows of the design matrix to the
        possible alternatives for this dataset.
    rows_to_mixers : 2D scipy sparse array.
        All elements should be zeros and ones. Will map the rows of the design
        matrix to the particular units that the mixing is being performed over.
        Note that in the case of panel data, this matrix will be different from
        `rows_to_obs`.
    choice_vector : 1D ndarray.
        All elements should be either ones or zeros. There should be one row
        per observation per available alternative for the given observation.
        Elements denote the alternative which is chosen by the given
        observation with a 1 and a zero otherwise.
    utility_transform :  callable.
        Should accept a 1D array of systematic utility values, a 1D array of
        alternative IDs, and miscellaneous args and kwargs. Should return a 2D
        array whose elements contain the appropriately transformed systematic
        utility values, based on the current model being evaluated and the
        given draw of the random coefficients. There should be one column for
        each draw of the random coefficients. There should have one row per
        individual per choice situation per available alternative.
    ridge : int, float, long, or None, optional.
        Determines whether or not ridge regression is performed. If a float is
        passed, then that float determines the ridge penalty for the
        optimization. Default = None.
    weights : 1D ndarray or None, optional.
        Allows for the calculation of weighted log-likelihoods. The weights can
        represent various things. In stratified samples, the weights may be
        the proportion of the observations in a given strata for a sample in
        relation to the proportion of observations in that strata in the
        population. In latent class models, the weights may be the probability
        of being a particular class. Default == None.

    Returns
    -------
    bhhh_matrix : 2D ndarray of shape `(design.shape[1], design.shape[1])`.
        The returned array is the BHHH approximation of the Fisher Information
        Matrix. I.e it is the negative of the sum of the outer product of
        each individual's gradient with itself.