docs: def calc_mixed_logit_gradient(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.
          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.
  
      Returns
      -------
      gradient : ndarray of shape (design_3d.shape[2],).
          The returned array is the gradient of the log-likelihood of the mixed
          MNL model with respect to `params`.
      """
      # Calculate the weights for the sample
      if weights is None:
          weights = np.ones(design_3d.shape[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 = (scaled_error[:, :, None] *
                  design_3d *
                  weights[:, None, None]).sum(axis=0)
  
      gradient = gradient.mean(axis=0)
  
      # Account for the ridge parameter if an L2 penalization is being performed
      if ridge is not None:
          gradient -= 2 * ridge * params
  
      return gradient.ravel()